R Markdown: The Definitive Guide https://bookdown.org/yihui/rmarkdown/

Tools for Reproducible Research http://kbroman.org/Tools4RR/

options(scipen = 999)  ## disable scientific notation 4.835312e-04 --> 0.0004835312

library(tidyverse) 
#> ✔ ggplot2 2.2.1.9000     ✔ purrr   0.2.4     
#> ✔ tibble  1.4.2          ✔ dplyr   0.7.4     
#> ✔ tidyr   0.8.0          ✔ stringr 1.3.0     
#> ✔ readr   1.1.1          ✔ forcats 0.3.0
library(readxl) #read excel
library(stargazer) #nice and informative tables
library(compareGroups) # nice table one
library(summarytools) # summary table for exploratory analysis
library(labelVector) # labeling variables
library(PoEdata)#for PoE datasets

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot. You can use include=FALSE to exclude everything in a chunk. If you only want to suppress messages, use message=FALSE instead. We can mask the result by using results = FALSE. Use command-shift-m to insert pipe operator %>%.

0.1 Useful tricks

#rename several variables
#dplyr::rename(.data, NEW.NAME1 = old.name1, NEW.NAME2 = old.name2)


## An 'if' statement takes the form:
## if(logical boolean){function to perform if the boolean is TRUE}
x <- 100000000

if(x > 10000){
  "x is a really big number"
}
## [1] "x is a really big number"
if(x < 10000){
  "x isn't that big"
} ## Nothing happens here because x < 10000 evaluates to false

## We can use else statements to deal with instances where the boolean
## evaluates to FALSE

## Note that using multiple if/else statements requires that you wrap the whole thing
## in curly brackets, {}:

x <- 1

if(x > 10000){
  "x is a really big number"
  }  else{
    "x isn't that big"
  }
## [1] "x isn't that big"
#######################################################

## Note about the modulus function: what it does is divides the number
## in front of it by the number after it, and gives you the remainder.

## It's really useful to use inside of loops because you can use it to
## print the iteration you're on.

## In other words, if you have a loop that has 100,000 iterations, 
## maybe you want to keep tabs on how long it's taking and where it is,
## but don't want to print every single number from 1 to 100,000.
## Inserting this code into the loop will print only every 100th iteration:

for(i in 1:500){
  if((i %% 100) == 0){
    print(i)
  }
}
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
####################################################

## For a loop that does 10,000 iterations:
start <- proc.time()[3] ## this captures and stores the time before the loop starts
holder <- c()

for(i in 1:10000){
  holder[i] <- rnorm(1, mean = 0, sd = 1) 
  # here we are drawing one observation from a normal distribution
  # with mean zero and standard deviation 1
}

end <- proc.time()[3]

## How long did that take?

loop10000 <- end - start
loop10000
## elapsed 
##   0.035
## What if we had vectorized our code and drawn all 10,000 random numbers first?

start <- proc.time()[3]
draws <-  rnorm(10000, mean = 0, sd = 1)
end <- proc.time()[3]

vect10000 <- end - start

loop10000
## elapsed 
##   0.035
vect10000
## elapsed 
##   0.001
## The vectorized code was way, way faster!

####################################################

## There's one other way, and it's a really cool R-hack.
## When you write a function, you can store objects using
## <<- instead of just <-
## Using the double arrow stores the object to your workspace,
## so that it's accessible after the function is complete

debugging.example <- function(x){
  
  holder <- rbeta(1000, 1, 10)
  holder <- sort(holder, decreasing = T)
  holder <- matrix(holder[1:50], 25, 2)
  # put the <<- here to store holder
  holder <<- matrix(holder[1:50], 25, 2)
  t(holder) %*% x
  
}

#debugging.example(test.matrix) ## Now it works
head(holder)
## [1]  0.38303320 -0.49538108 -0.14127816  0.01735916  0.64212009 -0.96262673
# Another trick to debug functions is to store objects in it with
# '<<-' instead of '<-'. If you only use '<-' then the only thing
# left in your workspace after the function runs will be the output
# from the final line of the function. Using '<<-' will force that
# object into your workspace, even if it's not the last line in
# the function

####################################################
# Note that with lists you can use lapply() or sapply()
# lapply() will output your result as a list
# sapply() will (usually) output your result as a vector
####################################################
## Debugging functions

## The browser() function is super super helpful if you have a really
## complicated function that has a bug in the code.

## Now let's see how we could use the browser() function to debug.
## Note: you can also put browser() functions into for/while loops


debugging.example <- function(x){
  
  holder <- rbeta(1000, 1, 10)
  browser()
  
  holder <- sort(holder, decreasing = T)
  browser()
  
  holder <- matrix(holder[1:50], 25, 2)
  browser()
  
  holder %*% x
}

test.matrix <- matrix(rnorm(100), 25, 4)
debugging.example(test.matrix) 
## this tells us that the error doesn't occur until after the 3rd browser() statement

0.2 Parallel computing

library(parallel)
library(doParallel)
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## Loading required package: iterators
library(foreach)

#Starting a cluster
numCores <- detectCores()
numCores
## [1] 8
cl <- makeCluster(numCores-1)





#Close the cluster
stopCluster(cl)

0.3 Exploratory analysis

#Introduction to summarytools
#https://cran.r-project.org/web/packages/summarytools/vignettes/Introduction.html

#install.packages('devtools')
#library(devtools)
#install_github('dcomtois/summarytools')

#view(dfSummary(iris))    #open a new browser to view the summary
#graph.magnif-> graph size
print(dfSummary(tobacco, style = 'grid', plain.ascii = FALSE, graph.magnif = 4), 
       method = 'render', headings = FALSE)
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 gender [factor] 1. F 2. M 489 (50.0%) 489 (50.0%) 978 (97.8%) 22 (2.2%)
2 age [numeric] mean (sd) : 49.6 (18.29) min < med < max : 18 < 50 < 80 IQR (CV) : 32 (0.37) 63 distinct values 975 (97.5%) 25 (2.5%)
3 age.gr [factor] 1. 18-34 2. 35-50 3. 51-70 4. 71 + 258 (26.5%) 241 (24.7%) 317 (32.5%) 159 (16.3%) 975 (97.5%) 25 (2.5%)
4 BMI [numeric] mean (sd) : 25.73 (4.49) min < med < max : 8.83 < 25.62 < 39.44 IQR (CV) : 5.72 (0.17) 974 distinct values 974 (97.4%) 26 (2.6%)
5 smoker [factor] 1. Yes 2. No 298 (29.8%) 702 (70.2%) 1000 (100%) 0 (0%)
6 cigs.per.day [numeric] mean (sd) : 6.78 (11.88) min < med < max : 0 < 0 < 40 IQR (CV) : 11 (1.75) 37 distinct values 965 (96.5%) 35 (3.5%)
7 diseased [factor] 1. Yes 2. No 224 (22.4%) 776 (77.6%) 1000 (100%) 0 (0%)
8 disease [character] 1. Hypertension 2. Cancer 3. Cholesterol 4. Heart 5. Pulmonary 6. Musculoskeletal 7. Diabetes 8. Hearing 9. Digestive 10. Hypotension [ 3 others ] 36 (16.2%) 34 (15.3%) 21 (9.5%) 20 (9.0%) 20 (9.0%) 19 (8.6%) 14 (6.3%) 14 (6.3%) 12 (5.4%) 11 (5.0%) 21 (9.5%) 222 (22.2%) 778 (77.8%)
9 samp.wgts [numeric] mean (sd) : 1 (0.08) min < med < max : 0.86 < 1.04 < 1.06 IQR (CV) : 0.19 (0.08) 0.86! : 267 (26.7%) 1.04! : 249 (24.9%) 1.05! : 324 (32.4%) 1.06! : 160 (16.0%) ! rounded 1000 (100%) 0 (0%)

Generated by summarytools 0.8.9 (R version 3.5.2)
2019-01-07

#http://lutgw1.lunet.edu/usr/lib64/R/library/Hmisc/tests/examples.Rmd
#http://hbiostat.org/R/Hmisc/examples.html
library(Hmisc)
d<-describe(tobacco) 
html(d, size=90, scroll=F)
tobacco

9 Variables   1000 Observations
gender
nmissingdistinct
978222
 Value        F   M
 Frequency  489 489
 Proportion 0.5 0.5
age
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
9752563149.621.1221.024.034.050.066.074.077.3
lowest : 18 19 20 21 22 , highest: 76 77 78 79 80
age.gr
image
nmissingdistinct
975254
 Value      18-34 35-50 51-70  71 +
 Frequency    258   241   317   159
 Proportion 0.265 0.247 0.325 0.163
BMI
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
97426974125.735.03918.4720.2422.9325.6228.6531.5533.11
lowest : 8.825777 9.00871310.34613112.98415714.179153
highest:38.12207438.37224438.63769439.20775339.439012
smoker
nmissingdistinct
100002
 Value        Yes    No
 Frequency    298   702
 Proportion 0.298 0.702
cigs.per.day
image
nmissingdistinctInfoMeanGmd.05.10.25.50.75.90.95
96535370.6596.78210.6 0 0 0 0112935
lowest : 0 5 6 7 8 , highest: 36 37 38 39 40
diseased
nmissingdistinct
100002
 Value        Yes    No
 Frequency    224   776
 Proportion 0.224 0.776
disease
image
nmissingdistinct
22277813
Cancer (34, 0.153), Cholesterol (21, 0.095), Diabetes (14, 0.063), Digestive (12, 0.054), Hearing (14, 0.063), Heart (20, 0.090), Hypertension (36, 0.162), Hypotension (11, 0.050), Musculoskeletal (19, 0.086), Neurological (10, 0.045), Other (2, 0.009), Pulmonary (20, 0.090), Vision (9, 0.041)
samp.wgts
image
nmissingdistinctInfoMeanGmd
1000040.92710.07774
 Value      0.8614232 1.0441767 1.0493827 1.0625000
 Frequency        267       249       324       160
 Proportion     0.267     0.249     0.324     0.160

0.4 Descriptive table one

#compareGroups 4.0: Descriptives by groups
#https://cran.r-project.org/web/packages/compareGroups/vignettes/compareGroups_vignette.html

data(predimed)
predimed$tmain <- with(predimed, Surv(toevent, event == "Yes"))
Hmisc::label(predimed$tmain) <- "AMI, stroke, or CV Death"

# https://cran.r-project.org/web/packages/labelVector/vignettes/labelVector.html
labelVector::get_label(predimed)
##  [1] "Intervention group"              "Sex"                            
##  [3] "Age"                             "Smoking"                        
##  [5] "Body mass index"                 "Waist circumference"            
##  [7] "Waist-to-height ratio"           "Hypertension"                   
##  [9] "Type-2 diabetes"                 "Dyslipidemia"                   
## [11] "Family history of premature CHD" "Hormone-replacement therapy"    
## [13] "MeDiet Adherence score"          "follow-up to main event (years)"
## [15] "AMI, stroke, or CV Death"        "AMI, stroke, or CV Death"
#set label
x <- 1:10
x <- set_label(x, "some integers")
x
## some integers
##  [1]  1  2  3  4  5  6  7  8  9 10
#To remove the variable toevent and event from the analysis:
compareGroups(group ~ . - toevent - event, data = predimed)
## 
## 
## -------- Summary of results by groups of 'Intervention group'---------
## 
## 
##    var                             N    p.value  method           
## 1  Sex                             6324 <0.001** categorical      
## 2  Age                             6324 0.003**  continuous normal
## 3  Smoking                         6324 0.444    categorical      
## 4  Body mass index                 6324 <0.001** continuous normal
## 5  Waist circumference             6324 0.045**  continuous normal
## 6  Waist-to-height ratio           6324 <0.001** continuous normal
## 7  Hypertension                    6324 0.249    categorical      
## 8  Type-2 diabetes                 6324 0.017**  categorical      
## 9  Dyslipidemia                    6324 0.423    categorical      
## 10 Family history of premature CHD 6324 0.581    categorical      
## 11 Hormone-replacement therapy     5661 0.850    categorical      
## 12 MeDiet Adherence score          6324 <0.001** continuous normal
## 13 AMI, stroke, or CV Death        6324 0.011**  Surv [Tmax=4.79] 
##    selection
## 1  ALL      
## 2  ALL      
## 3  ALL      
## 4  ALL      
## 5  ALL      
## 6  ALL      
## 7  ALL      
## 8  ALL      
## 9  ALL      
## 10 ALL      
## 11 ALL      
## 12 ALL      
## 13 ALL      
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1

No filters have been used (e.g., selecting only treated patients); therefore, the selection column lists ALL (for all variables).

#To perform the analysis in a subset of participants (e.g., “female” participants):
compareGroups(group ~ age + smoke + waist + hormo, data = predimed, 
    subset = sex == "Female")
## 
## 
## -------- Summary of results by groups of 'group'---------
## 
## 
##   var                         N    p.value method           
## 1 Age                         3645 0.056*  continuous normal
## 2 Smoking                     3645 0.907   categorical      
## 3 Waist circumference         3645 0.016** continuous normal
## 4 Hormone-replacement therapy 3459 0.898   categorical      
##   selection      
## 1 sex == "Female"
## 2 sex == "Female"
## 3 sex == "Female"
## 4 sex == "Female"
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#To subset specific variable/s (e.g., hormo and waist):
compareGroups(group ~ age + sex + smoke + waist + hormo, data = predimed, 
    selec = list(hormo = sex == "Female", waist = waist > 20))
## 
## 
## -------- Summary of results by groups of 'Intervention group'---------
## 
## 
##   var                         N    p.value  method           
## 1 Age                         6324 0.003**  continuous normal
## 2 Sex                         6324 <0.001** categorical      
## 3 Smoking                     6324 0.444    categorical      
## 4 Waist circumference         6324 0.045**  continuous normal
## 5 Hormone-replacement therapy 3459 0.898    categorical      
##   selection      
## 1 ALL            
## 2 ALL            
## 3 ALL            
## 4 waist > 20     
## 5 sex == "Female"
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#Combinations are also allowed, e.g.:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed, 
    selec = list(waist = !is.na(hormo)), subset = sex == "Female")
## 
## 
## -------- Summary of results by groups of 'group'---------
## 
## 
##   var                         N    p.value method           
## 1 Age                         3645 0.056*  continuous normal
## 2 Smoking                     3645 0.907   categorical      
## 3 Waist circumference         3459 0.007** continuous normal
## 4 Hormone-replacement therapy 3459 0.898   categorical      
##   selection                          
## 1 sex == "Female"                    
## 2 sex == "Female"                    
## 3 (sex == "Female") & (!is.na(hormo))
## 4 sex == "Female"                    
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#A variable can appear twice in the formula, e.g.:  (NOT SHOW)

#To change default options, e.g., “waist” used as non-normal distributed:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed, 
    method = c(waist = 2))
## Warning in cor.test.default(x, as.integer(y), method = "spearman"): Cannot
## compute exact p-value with ties
## 
## 
## -------- Summary of results by groups of 'Intervention group'---------
## 
## 
##   var                         N    p.value method                selection
## 1 Age                         6324 0.003** continuous normal     ALL      
## 2 Smoking                     6324 0.444   categorical           ALL      
## 3 Waist circumference         6324 0.085*  continuous non-normal ALL      
## 4 Hormone-replacement therapy 5661 0.850   categorical           ALL      
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#If the method argument is stated as NA for a variable, then a Shapiro-Wilk test for normality is used to decide if the variable is normal or non-normal distributed. To change the significance threshold:
compareGroups(group ~ age + smoke + waist + hormo, data = predimed, 
    method = c(waist = NA), alpha = 0.01)
## Warning in cor.test.default(x, as.integer(y), method = "spearman"): Cannot
## compute exact p-value with ties
## 
## 
## -------- Summary of results by groups of 'Intervention group'---------
## 
## 
##   var                         N    p.value method                selection
## 1 Age                         6324 0.003** continuous normal     ALL      
## 2 Smoking                     6324 0.444   categorical           ALL      
## 3 Waist circumference         6324 0.085*  continuous non-normal ALL      
## 4 Hormone-replacement therapy 5661 0.850   categorical           ALL      
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#Shapiro-Wilk test, stating the cutpoint at 0.01 level
#1: forces analysis as normal-distributed
#2: forces analysis as continuous non-normal
#3: forces analysis as categorical
#NA: performs a Shapiro-Wilks test to decide between normal or non-normal

#All non factor variables are considered as continuous. Exception is made (by default) for those that have fewer than 5 different values. 

#The object from compareGroups can later be updated. For example:
res <- compareGroups(group ~ age + sex + smoke + waist + hormo, 
    data = predimed)
res <- update(res, . ~ . - sex + bmi + toevent, subset = sex == 
    "Female", method = c(waist = 2, tovent = 2), selec = list(bmi = !is.na(hormo)))
## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## variables tovent specified in 'method' not found

## Warning in compareGroups.fit(X = X, y = y, include.label = include.label, :
## Cannot compute exact p-value with ties
res
## 
## 
## -------- Summary of results by groups of 'group'---------
## 
## 
##   var                             N    p.value  method               
## 1 Age                             3645 0.056*   continuous normal    
## 2 Smoking                         3645 0.907    categorical          
## 3 Waist circumference             3645 0.037**  continuous non-normal
## 4 Hormone-replacement therapy     3459 0.898    categorical          
## 5 Body mass index                 3459 0.002**  continuous normal    
## 6 follow-up to main event (years) 3645 <0.001** continuous normal    
##   selection                          
## 1 sex == "Female"                    
## 2 sex == "Female"                    
## 3 sex == "Female"                    
## 4 sex == "Female"                    
## 5 (sex == "Female") & (!is.na(hormo))
## 6 sex == "Female"                    
## -----
## Signif. codes:  0 '**' 0.05 '*' 0.1 ' ' 1
#the Odds Ratio (OR) can be printed in the final table. If the response variable is time-to-event (see Section 3.1), the Hazard Ratio (HR) can be printed instead.
res2 <- compareGroups(htn ~ age + sex + bmi + smoke, data = predimed, 
    ref = c(smoke = 1, sex = 2))
createTable(res2, show.ratio = TRUE)
## 
## --------Summary descriptives table by 'Hypertension'---------
## 
## ___________________________________________________________________________ 
##                     No          Yes             OR        p.ratio p.overall 
##                   N=1089       N=5235                                       
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Age             65.9 (6.19) 67.2 (6.15)  1.04 [1.03;1.05] <0.001   <0.001   
## Sex:                                                               <0.001   
##     Male        595 (54.6%) 2084 (39.8%) 0.55 [0.48;0.63]  0.000            
##     Female      494 (45.4%) 3151 (60.2%)       Ref.        Ref.             
## Body mass index 28.9 (3.69) 30.2 (3.80)  1.10 [1.08;1.12] <0.001   <0.001   
## Smoking:                                                           <0.001   
##     Never       536 (49.2%) 3356 (64.1%)       Ref.        Ref.             
##     Current     233 (21.4%) 625 (11.9%)  0.43 [0.36;0.51]  0.000            
##     Former      320 (29.4%) 1254 (24.0%) 0.63 [0.54;0.73] <0.001            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#The createTable function, applied to an object of compareGroups class, returns tables with descriptives that can be displayed on-screen or exported to CSV, LaTeX, HTML, Word or Excel.
res <- compareGroups(group ~ age + sex + smoke + waist + hormo, 
    data = predimed, selec = list(hormo = sex == "Female"))
restab <- createTable(res)
#the option “descr” prints descriptive tables.
#digits: The number of digits that appear in the results can be changed
print(restab, which.table = "descr")
## 
## --------Summary descriptives table by 'Intervention group'---------
## 
## ________________________________________________________________________________ 
##                                Control    MedDiet + Nuts MedDiet + VOO p.overall 
##                                 N=2042        N=2100        N=2182               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Age                          67.3 (6.28)   66.7 (6.02)    67.0 (6.21)    0.003   
## Sex:                                                                    <0.001   
##     Male                     812 (39.8%)   968 (46.1%)    899 (41.2%)            
##     Female                   1230 (60.2%)  1132 (53.9%)  1283 (58.8%)            
## Smoking:                                                                 0.444   
##     Never                    1282 (62.8%)  1259 (60.0%)  1351 (61.9%)            
##     Current                  270 (13.2%)   296 (14.1%)    292 (13.4%)            
##     Former                   490 (24.0%)   545 (26.0%)    539 (24.7%)            
## Waist circumference           101 (10.8)    100 (10.6)    100 (10.4)     0.045   
## Hormone-replacement therapy:                                             0.898   
##     No                       1143 (97.4%)  1036 (97.2%)  1183 (97.0%)            
##     Yes                       31 (2.64%)    30 (2.81%)    36 (2.95%)             
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#show.all: If show.all argument is set to TRUE a column is displayed with descriptives for all data:
createTable(res, show.all = TRUE)
## 
## --------Summary descriptives table by 'Intervention group'---------
## 
## _____________________________________________________________________________________________ 
##                                 [ALL]       Control    MedDiet + Nuts MedDiet + VOO p.overall 
##                                 N=6324       N=2042        N=2100        N=2182               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Age                          67.0 (6.17)  67.3 (6.28)   66.7 (6.02)    67.0 (6.21)    0.003   
## Sex:                                                                                 <0.001   
##     Male                     2679 (42.4%) 812 (39.8%)   968 (46.1%)    899 (41.2%)            
##     Female                   3645 (57.6%) 1230 (60.2%)  1132 (53.9%)  1283 (58.8%)            
## Smoking:                                                                              0.444   
##     Never                    3892 (61.5%) 1282 (62.8%)  1259 (60.0%)  1351 (61.9%)            
##     Current                  858 (13.6%)  270 (13.2%)   296 (14.1%)    292 (13.4%)            
##     Former                   1574 (24.9%) 490 (24.0%)   545 (26.0%)    539 (24.7%)            
## Waist circumference           100 (10.6)   101 (10.8)    100 (10.6)    100 (10.4)     0.045   
## Hormone-replacement therapy:                                                          0.898   
##     No                       3362 (97.2%) 1143 (97.4%)  1036 (97.2%)  1183 (97.0%)            
##     Yes                       97 (2.80%)   31 (2.64%)    30 (2.81%)    36 (2.95%)             
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Tables made with the same response variable can be combined by row:
restab1 <- createTable(compareGroups(group ~ age + sex, data = predimed))
restab2 <- createTable(compareGroups(group ~ bmi + smoke, data = predimed))
rbind(`Non-modifiable risk factors` = restab1, `Modifiable risk factors` = restab2)
## 
## --------Summary descriptives table by 'Intervention group'---------
## 
## _______________________________________________________________________ 
##                       Control    MedDiet + Nuts MedDiet + VOO p.overall 
##                        N=2042        N=2100        N=2182               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Non-modifiable risk factors:
##     Age             67.3 (6.28)   66.7 (6.02)    67.0 (6.21)    0.003   
##     Sex:                                                       <0.001   
##         Male        812 (39.8%)   968 (46.1%)    899 (41.2%)            
##         Female      1230 (60.2%)  1132 (53.9%)  1283 (58.8%)            
## Modifiable risk factors:
##     Body mass index 30.3 (3.96)   29.7 (3.77)    29.9 (3.71)   <0.001   
##     Smoking:                                                    0.444   
##         Never       1282 (62.8%)  1259 (60.0%)  1351 (61.9%)            
##         Current     270 (13.2%)   296 (14.1%)    292 (13.4%)            
##         Former      490 (24.0%)   545 (26.0%)    539 (24.7%)            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Columns from tables built with the same explanatory and response variables but done with a different subset (i.e. “Male” and “Female”, strata) can be combined:

# first build the table with descriptives by groups:
# include.miss = T shows missing value category 
res <- compareGroups(group ~ . - sex, include.miss = T, predimed)
restab <- createTable(res, hide.no = "no")
# and then apply the strataTable function on the table:
strataTable(restab, "sex")
## 
## --------Summary descriptives table ---------
## 
## _______________________________________________________________________________________________________________________________________
##                                                        Male                                               Female                       
##                                 __________________________________________________  ___________________________________________________
##                                   Control   MedDiet + Nuts MedDiet + VOO p.overall    Control    MedDiet + Nuts MedDiet + VOO p.overall 
##                                    N=812        N=968          N=899                   N=1230        N=1132        N=1283               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## Age                             66.4 (6.62)  65.8 (6.40)    66.1 (6.61)    0.215    68.0 (5.96)   67.4 (5.57)    67.7 (5.84)    0.056   
## Smoking:                                                                   0.851                                                0.907   
##     Never                       205 (25.2%)  266 (27.5%)    236 (26.3%)             1077 (87.6%)  993 (87.7%)   1115 (86.9%)            
##     Current                     204 (25.1%)  242 (25.0%)    221 (24.6%)              66 (5.37%)    54 (4.77%)    71 (5.53%)             
##     Former                      403 (49.6%)  460 (47.5%)    442 (49.2%)              87 (7.07%)    85 (7.51%)    97 (7.56%)             
## Body mass index                 29.6 (3.45)  29.1 (3.28)    29.2 (3.28)    0.018    30.8 (4.20)   30.2 (4.08)    30.4 (3.91)    0.002   
## Waist circumference             104 (9.82)    103 (9.36)    103 (9.65)     0.289    99.0 (11.0)   97.8 (11.0)    98.0 (10.5)    0.016   
## Waist-to-height ratio           0.62 (0.06)  0.62 (0.06)    0.62 (0.06)    0.191    0.64 (0.07)   0.63 (0.07)    0.63 (0.07)    0.002   
## Hypertension                    649 (79.9%)  753 (77.8%)    682 (75.9%)    0.130    1062 (86.3%)  985 (87.0%)   1104 (86.0%)    0.780   
## Type-2 diabetes                 430 (53.0%)  496 (51.2%)    486 (54.1%)    0.468    540 (43.9%)   454 (40.1%)    596 (46.5%)    0.007   
## Dyslipidemia                    531 (65.4%)  653 (67.5%)    589 (65.5%)    0.575    948 (77.1%)   886 (78.3%)    971 (75.7%)    0.319   
## Family history of premature CHD 135 (16.6%)  171 (17.7%)    156 (17.4%)    0.841    327 (26.6%)   289 (25.5%)    351 (27.4%)    0.596   
## Hormone-replacement therapy:                                               0.903                                                0.689   
##     No                          668 (82.3%)  799 (82.5%)    735 (81.8%)             1143 (92.9%)  1036 (91.5%)  1183 (92.2%)            
##     Yes                          0 (0.00%)    0 (0.00%)      0 (0.00%)               31 (2.52%)    30 (2.65%)    36 (2.81%)             
##     'Missing'                   144 (17.7%)  169 (17.5%)    164 (18.2%)              56 (4.55%)    66 (5.83%)    64 (4.99%)             
## MeDiet Adherence score          8.57 (1.94)  8.86 (1.96)    8.91 (2.02)    0.001    8.36 (1.94)   8.77 (1.86)    8.68 (1.92)   <0.001   
## follow-up to main event (years) 4.05 (1.78)  4.38 (1.74)    4.53 (1.64)   <0.001    4.12 (1.71)   4.26 (1.67)    4.72 (1.56)   <0.001   
## AMI, stroke, or CV Death        58 (7.14%)    41 (4.24%)    52 (5.78%)     0.029     39 (3.17%)    29 (2.56%)    33 (2.57%)     0.576   
## AMI, stroke, or CV Death           8.48%        4.42%          5.71%       0.008       3.84%         2.66%          2.16%       0.269   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
export2md(strataTable(restab, "sex"), size = 8)
Summary descriptive tables

Male
Female
Control MedDiet + Nuts MedDiet + VOO p.overall Control MedDiet + Nuts MedDiet + VOO p.overall
N=812 N=968 N=899 N=1230 N=1132 N=1283
Age 66.4 (6.62) 65.8 (6.40) 66.1 (6.61) 0.215 68.0 (5.96) 67.4 (5.57) 67.7 (5.84) 0.056
Smoking: 0.851 0.907
Never 205 (25.2%) 266 (27.5%) 236 (26.3%) 1077 (87.6%) 993 (87.7%) 1115 (86.9%)
Current 204 (25.1%) 242 (25.0%) 221 (24.6%) 66 (5.37%) 54 (4.77%) 71 (5.53%)
Former 403 (49.6%) 460 (47.5%) 442 (49.2%) 87 (7.07%) 85 (7.51%) 97 (7.56%)
Body mass index 29.6 (3.45) 29.1 (3.28) 29.2 (3.28) 0.018 30.8 (4.20) 30.2 (4.08) 30.4 (3.91) 0.002
Waist circumference 104 (9.82) 103 (9.36) 103 (9.65) 0.289 99.0 (11.0) 97.8 (11.0) 98.0 (10.5) 0.016
Waist-to-height ratio 0.62 (0.06) 0.62 (0.06) 0.62 (0.06) 0.191 0.64 (0.07) 0.63 (0.07) 0.63 (0.07) 0.002
Hypertension 649 (79.9%) 753 (77.8%) 682 (75.9%) 0.130 1062 (86.3%) 985 (87.0%) 1104 (86.0%) 0.780
Type-2 diabetes 430 (53.0%) 496 (51.2%) 486 (54.1%) 0.468 540 (43.9%) 454 (40.1%) 596 (46.5%) 0.007
Dyslipidemia 531 (65.4%) 653 (67.5%) 589 (65.5%) 0.575 948 (77.1%) 886 (78.3%) 971 (75.7%) 0.319
Family history of premature CHD 135 (16.6%) 171 (17.7%) 156 (17.4%) 0.841 327 (26.6%) 289 (25.5%) 351 (27.4%) 0.596
Hormone-replacement therapy: 0.903 0.689
No 668 (82.3%) 799 (82.5%) 735 (81.8%) 1143 (92.9%) 1036 (91.5%) 1183 (92.2%)
Yes 0 (0.00%) 0 (0.00%) 0 (0.00%) 31 (2.52%) 30 (2.65%) 36 (2.81%)
‘Missing’ 144 (17.7%) 169 (17.5%) 164 (18.2%) 56 (4.55%) 66 (5.83%) 64 (4.99%)
MeDiet Adherence score 8.57 (1.94) 8.86 (1.96) 8.91 (2.02) 0.001 8.36 (1.94) 8.77 (1.86) 8.68 (1.92) <0.001
follow-up to main event (years) 4.05 (1.78) 4.38 (1.74) 4.53 (1.64) <0.001 4.12 (1.71) 4.26 (1.67) 4.72 (1.56) <0.001
AMI, stroke, or CV Death 58 (7.14%) 41 (4.24%) 52 (5.78%) 0.029 39 (3.17%) 29 (2.56%) 33 (2.57%) 0.576
AMI, stroke, or CV Death 8.48% 4.42% 5.71% 0.008 3.84% 2.66% 2.16% 0.269 
#Tables can be exported to CSV, HTML, LaTeX, PDF, Markdown, Word or Excel
#export2csv(restab, file='table1.csv'), exports to CSV format
#export2html(restab, file='table1.html'), exports to HTML format
#export2latex(restab, file='table1.tex'), exports to LaTeX format (to be included in Swaeave documents R chunks)
#export2pdf(restab, file='table1.pdf'), exports to PDF format
#export2md(restab, file='table1.md'), to be included inside Markdown documents R chunks
#export2word(restab, file='table1.docx'), exports to Word format
#export2xls(restab, file='table1.xlsx'), exports to Excel format

0.5 Reporting and handling missing data

A new function has been implemented in the compareGroups package, which is called missingTable. This function applies to both compareGroups and createTable class objects. This last option is useful when the table is already created. To illustrate it, we will use the REGICOR data set, comparing missing rates of all variables by year:

# from a compareGroups object
data(regicor)
res <- compareGroups(year ~ . - id, regicor)
missingTable(res)
## Warning in cor(as.integer(x), as.integer(y)): the standard deviation is
## zero

## Warning in cor(as.integer(x), as.integer(y)): the standard deviation is
## zero
## 
## --------Missingness table by 'year'---------
## 
## ____________________________________________________ 
##             1995       2000        2005    p.overall 
##            N=431       N=786      N=1077             
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## age      0 (0.00%)   0 (0.00%)  0 (0.00%)      .     
## sex      0 (0.00%)   0 (0.00%)  0 (0.00%)      .     
## smoker   16 (3.71%) 28 (3.56%)  17 (1.58%)   0.010   
## sbp      3 (0.70%)  11 (1.40%)  0 (0.00%)   <0.001   
## dbp      3 (0.70%)  11 (1.40%)  0 (0.00%)   <0.001   
## histhtn  0 (0.00%)   0 (0.00%)  8 (0.74%)    0.015   
## txhtn    0 (0.00%)   0 (0.00%)  43 (3.99%)  <0.001   
## chol     28 (6.50%) 71 (9.03%)  2 (0.19%)   <0.001   
## hdl      30 (6.96%) 38 (4.83%)  1 (0.09%)   <0.001   
## triglyc  28 (6.50%) 34 (4.33%)  1 (0.09%)   <0.001   
## ldl      43 (9.98%) 98 (12.5%)  27 (2.51%)  <0.001   
## histchol 0 (0.00%)  15 (1.91%)  6 (0.56%)    0.001   
## txchol   0 (0.00%)  13 (1.65%)  42 (3.90%)  <0.001   
## height   8 (1.86%)  15 (1.91%)  12 (1.11%)   0.318   
## weight   8 (1.86%)  15 (1.91%)  12 (1.11%)   0.318   
## bmi      8 (1.86%)  15 (1.91%)  12 (1.11%)   0.318   
## phyact   64 (14.8%) 22 (2.80%)  2 (0.19%)   <0.001   
## pcs      34 (7.89%) 123 (15.6%) 83 (7.71%)  <0.001   
## mcs      34 (7.89%) 123 (15.6%) 83 (7.71%)  <0.001   
## cv       33 (7.66%) 45 (5.73%)  53 (4.92%)   0.118   
## tocv     33 (7.66%) 45 (5.73%)  53 (4.92%)   0.118   
## death    44 (10.2%) 48 (6.11%)  54 (5.01%)   0.001   
## todeath  44 (10.2%) 48 (6.11%)  54 (5.01%)   0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#summarize missingness
# the funcation counts the # of missing value for each variable
propmiss <- function(dataframe) {
m <- sapply(dataframe, function(x) {
data.frame(
nmiss=sum(is.na(x)), 
n=length(x), 
propmiss=sum(is.na(x))/length(x)
        )
    })
d <- data.frame(t(m))
d <- sapply(d, unlist)
d <- as.data.frame(d)
d$variable <- row.names(d)
    row.names(d) <- NULL
d <- cbind(d[ncol(d)],d[-ncol(d)])
return(d[order(d$propmiss), ])
}


#replace NA to mean and mode
Mode<-function(x,na.rm){
    xtab <- table(x)
    xmode <- names(which(xtab == max(xtab)))
    if (length(xmode) > 1) xmode <- ">1 mode"
    return(xmode)
}

df_test

for (var in 1:ncol(df_test)) {
    if (class(df_test[,var])=="numeric") {
       df_test[is.na(df_test[,var]),var] <- mean(df_test[,var], na.rm = TRUE)
    } else if (class(df_test[,var]) %in% c("character", "factor")) {
        df_test[is.na(df_test[,var]),var] <- Mode(df_test[,var], na.rm = TRUE)
    }
 }

ac<-df_test

0.6 Creating regression output table and export -> html -> excel

#load package with the data
#install.packages("devtools")  # if not already installed
# library(devtools)
# install_git("https://github.com/ccolonescu/PoEdata")

#load data 
data("cps4_small", package="PoEdata")
names(cps4_small)
##  [1] "wage"    "educ"    "exper"   "hrswk"   "married" "female"  "metro"  
##  [8] "midwest" "south"   "west"    "black"   "asian"
#run the model
dnosouth <- cps4_small[which(cps4_small$south==0),]#no south
dsouth <- cps4_small[which(cps4_small$south==1),] #south
mod5ns <- lm(wage~educ+black*female, data=dnosouth)
mod5s <- lm(wage~educ+black*female, data=dsouth)
mod6 <- lm(wage~educ+black*female+south/(educ+black*female),
                       data=cps4_small)

#Stargazer cheat sheets
#https://www.jakeruss.com/cheatsheets/stargazer/

print(stargazer(mod6, mod5ns, mod5s, header=FALSE, 
  type="text",
  title="Model comparison, 'wage' equation",
  keep.stat="n",digits=2, single.row=TRUE,
  intercept.bottom=FALSE))
## 
## Model comparison, 'wage' equation
## ==================================================================
##                                  Dependent variable:              
##                    -----------------------------------------------
##                                         wage                      
##                          (1)             (2)             (3)      
## ------------------------------------------------------------------
## Constant           -6.61*** (2.34) -6.61*** (2.30)  -2.66 (3.42)  
## educ               2.17*** (0.17)  2.17*** (0.16)  1.86*** (0.24) 
## black               -5.09* (2.64)   -5.09* (2.60)   -3.38 (2.58)  
## female             -5.01*** (0.90) -5.01*** (0.89) -4.10*** (1.58)
## south                3.94 (4.05)                                  
## black:female         5.31 (3.50)     5.31 (3.45)     2.37 (3.38)  
## educ:south          -0.31 (0.29)                                  
## black:south          1.70 (3.63)                                  
## female:south         0.90 (1.77)                                  
## black:female:south  -2.94 (4.79)                                  
## ------------------------------------------------------------------
## Observations            1,000            704             296      
## ==================================================================
## Note:                                  *p<0.1; **p<0.05; ***p<0.01
##  [1] ""                                                                  
##  [2] "Model comparison, 'wage' equation"                                 
##  [3] "=================================================================="
##  [4] "                                 Dependent variable:              "
##  [5] "                   -----------------------------------------------"
##  [6] "                                        wage                      "
##  [7] "                         (1)             (2)             (3)      "
##  [8] "------------------------------------------------------------------"
##  [9] "Constant           -6.61*** (2.34) -6.61*** (2.30)  -2.66 (3.42)  "
## [10] "educ               2.17*** (0.17)  2.17*** (0.16)  1.86*** (0.24) "
## [11] "black               -5.09* (2.64)   -5.09* (2.60)   -3.38 (2.58)  "
## [12] "female             -5.01*** (0.90) -5.01*** (0.89) -4.10*** (1.58)"
## [13] "south                3.94 (4.05)                                  "
## [14] "black:female         5.31 (3.50)     5.31 (3.45)     2.37 (3.38)  "
## [15] "educ:south          -0.31 (0.29)                                  "
## [16] "black:south          1.70 (3.63)                                  "
## [17] "female:south         0.90 (1.77)                                  "
## [18] "black:female:south  -2.94 (4.79)                                  "
## [19] "------------------------------------------------------------------"
## [20] "Observations            1,000            704             296      "
## [21] "=================================================================="
## [22] "Note:                                  *p<0.1; **p<0.05; ***p<0.01"
#output the results to html, we create a HTML file that can be easily copied/pasted into Word/Excel.
#HTML output may look like meaningless code in the console. use results = F to mask them.
# add "'" in front of the cell contents in excel to avoid function
xxx<-stargazer(mod6, mod5ns, mod5s, header=FALSE, 
  type="html",
  title="Model comparison, 'wage' equation",
  keep.stat="n",digits=2, single.row=TRUE,
  intercept.bottom=FALSE, out= "xxx.html")

library(rio)
#export(as.data.frame(xxx), "xxx.csv")

0.7 JJHmisc package (TODO)

#library(devtools)
#install_github("johnjosephhorton/JJHmisc")
library(JJHmisc)
library(help = JJHmisc)
lsf.str("package:JJHmisc")
## AddTableNote : function (stargazer.object, out.file, note, adjust = 3, table.type = "table")  
## big.letter.theme : function (font.size)  
## clusterSE : function (cluster.name, model, df)  
## Codebook.Entry : function (var.name, label, note)  
## createLaTeXparameters : function (key.list, parameter.file.location, append = FALSE)  
## dzip : function (keys, values)  
## fix.sample.size.line : function (table.line)  
## flattenList : function (first.key, values)  
## genParamAdder : function (parameter.file)  
## get_field_names : function (con, tablename)  
## get.line : function (linenumber, file)  
## GetData : function (table.name, file.name, refresh)  
## HasChange : function (x, k)  
## insert.note : function (linenumber, note, file)  
## kill.lines : function (start, end, file)  
## latex.column.label : function (width, number, label, include.DV = TRUE)  
## latex.row.label : function (width, label, rule.pts = 0, new.lines = 0, new.line.char = FALSE)  
## NoteWrapper : function (x, text.width = 0.9)  
## pp.big.number : function (e)  
## regression.table : function (models, renames, out.file)  
## render : function (x, lst)  
## spanning.latex.column.label : function (num.columns, width, label)  
## starLabel : function (p)  
## tableContinuous2 : function (vars, weights = NA, subset = NA, group = NA, stats = c("n", 
##     "min", "q1", "median", "mean", "q3", "max", "s", "iqr", "na"), 
##     prec = 1, col.tit = NA, col.tit.font = c("bf", "", "sf", "it", 
##         "rm"), print.pval = c("none", "anova", "kruskal"), pval.bound = 10^-4, 
##     declare.zero = 10^-10, cap = "", lab = "", font.size = "footnotesize", 
##     longtable = TRUE, disp.cols = NA, nams = NA)  
## tableNominal2 : function (vars, weights = NA, subset = NA, group = NA, miss.cat = NA, 
##     print.pval = c("none", "fisher", "chi2"), pval.bound = 10^-4, 
##     fisher.B = 2000, vertical = TRUE, cap = "", lab = "", col.tit.font = c("bf", 
##         "", "sf", "it", "rm"), font.size = "footnotesize", longtable = TRUE, 
##     nams = NA, cumsum = TRUE)  
## wasItUsed : function (r.file, data.frame)  
## write.image : function (filename, g, width = 5, height = 5, format = "png")  
## writeImage : function (g, file.name.no.ext, width = 15, height = 8, path = "../../writeup/plots/", 
##     position = "h", line.width = "0.8", caption = "", notes = "", 
##     include.tex.wrapper = FALSE)
#read the source code of AddTableNote function
JJHmisc:::AddTableNote 
## function (stargazer.object, out.file, note, adjust = 3, table.type = "table") 
## {
##     "Adds a table note to a stargazer table. It assumes that the\n     last three lines of stargazer output are useless."
##     n <- length(stargazer.object)
##     write(stargazer.object[1:(n - adjust)], out.file)
##     write(c("\\end{tabular}"), out.file, append = TRUE)
##     write(note, out.file, append = TRUE)
##     write(c(paste0("\\end{", table.type, "}")), out.file, append = TRUE)
## }
## <bytecode: 0x7fef18b20928>
## <environment: namespace:JJHmisc>

0.8 list of some useful R functions

knitr::include_graphics("rFunctionsList.pdf")
Useful functions

Useful functions

https://tfse.mikewk.com

0.9 R and stan

0.9.1 Simple example

library(rstan)
library(bayesplot)
stancode <- 'data {real y_mean;} 
              parameters {real y;} 
              model {y ~ normal(y_mean,1);}'
mod <- stan_model(model_code = stancode)
simple.fit <- sampling(mod, data = list(y_mean = 0))
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 7e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.006417 seconds (Warm-up)
## Chain 1:                0.006959 seconds (Sampling)
## Chain 1:                0.013376 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.006231 seconds (Warm-up)
## Chain 2:                0.00617 seconds (Sampling)
## Chain 2:                0.012401 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.006364 seconds (Warm-up)
## Chain 3:                0.006119 seconds (Sampling)
## Chain 3:                0.012483 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 3e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.0061 seconds (Warm-up)
## Chain 4:                0.005758 seconds (Sampling)
## Chain 4:                0.011858 seconds (Total)
## Chain 4:
simple.fit2 <- sampling(mod, data = list(y_mean = 5))
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3e-06 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.00611 seconds (Warm-up)
## Chain 1:                0.005841 seconds (Sampling)
## Chain 1:                0.011951 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 3e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.006205 seconds (Warm-up)
## Chain 2:                0.006822 seconds (Sampling)
## Chain 2:                0.013027 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.006055 seconds (Warm-up)
## Chain 3:                0.005907 seconds (Sampling)
## Chain 3:                0.011962 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '5a0eac41591944c7d7b299fda75d726b' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 2e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.006037 seconds (Warm-up)
## Chain 4:                0.006337 seconds (Sampling)
## Chain 4:                0.012374 seconds (Total)
## Chain 4:

0.9.2 Write an stan model file and run the analysis

Comments are indicated by // in Stan. The write("model code", "file_name") bit allows us to write the Stan model in our R script and output the file to the working directory (or you can set a different file path).

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {
} // The posterior predictive distribution",

"stan_model1.stan")

# stanc("stan_model1.stan")

stan_model1 <- "stan_model1.stan"

#load data

# Adding stringsAsFactors = F means that numeric variables won't be
# read in as factors/categorical variables
seaice <- read.csv("seaice.csv", stringsAsFactors = F)
head(seaice)
##   year extent_north extent_south
## 1 1979       12.328       11.700
## 2 1980       12.337       11.230
## 3 1981       12.127       11.435
## 4 1982       12.447       11.640
## 5 1983       12.332       11.389
## 6 1984       11.910       11.454
#To explore the answer to that question, first we can make a figure.
plot(extent_north ~ year, pch = 20, data = seaice)

#preparing data for stan
x <- I(seaice$year - 1978) #rename the variables and index the years from 1 to 39. 
y <- seaice$extent_north
N <- length(seaice$year)
stan_data <- list(N = N, x = x, y = y)

#fit the model with stan
stan_model1.fit <- stan(file = stan_model1, 
            data = stan_data, 
            warmup = 500, iter = 1000, 
            chains = 4, cores = 2, thin = 1)

stan_model1.fit
## Inference for Stan model: stan_model1.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## alpha 12.55    0.00 0.08 12.40 12.50 12.55 12.61 12.69   849    1
## beta  -0.05    0.00 0.00 -0.06 -0.06 -0.05 -0.05 -0.05   958    1
## sigma  0.23    0.00 0.03  0.18  0.21  0.22  0.24  0.28   904    1
## lp__  37.48    0.04 1.19 34.32 36.91 37.80 38.36 38.87   797    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Jan  7 16:27:03 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
#look at the full posterior of our parameters by extracting them from the model object.
#extract() puts the posterior estimates for each parameter into a list.
posterior <- extract(stan_model1.fit)
str(posterior)
## List of 4
##  $ alpha: num [1:2000(1d)] 12.5 12.5 12.5 12.4 12.6 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ beta : num [1:2000(1d)] -0.0525 -0.0523 -0.0532 -0.05 -0.0563 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ sigma: num [1:2000(1d)] 0.245 0.231 0.246 0.215 0.275 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
##  $ lp__ : num [1:2000(1d)] 38.2 38.6 38.2 37.5 37 ...
##   ..- attr(*, "dimnames")=List of 1
##   .. ..$ iterations: NULL
#We can also look at the posterior densities & histograms.
stan_dens(stan_model1.fit)

stan_hist(stan_model1.fit)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#we can generate plots which indicate the mean parameter estimates and
#any credible intervals we may be interested in. 
#Note that the 95% credible intervals for the beta and 
#sigma parameters are very small, thus you only see the dots.
plot(stan_model1.fit, show_density = FALSE, ci_level = 0.5, 
     outer_level = 0.95, fill_color = "salmon")
## ci_level: 0.5 (50% intervals)
## outer_level: 0.95 (95% intervals)

#comparing with lm function

lm1 <- lm(y ~ x)
summary(lm1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.49925 -0.17713  0.04898  0.16923  0.32829 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 12.555888   0.071985   174.4 <0.0000000000000002 ***
## x           -0.054574   0.003137   -17.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2205 on 37 degrees of freedom
## Multiple R-squared:  0.8911, Adjusted R-squared:  0.8881 
## F-statistic: 302.7 on 1 and 37 DF,  p-value: < 0.00000000000000022
lm_alpha <- summary(lm1)$coeff[1]  # the intercept
lm_beta <- summary(lm1)$coeff[2]  # the slope
lm_sigma <- sigma(lm1)  # the residual error

plot(y ~ x, pch = 20)

abline(lm1, col = 2, lty = 2, lw = 3)
abline( mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)

The result is identical to the lm output. This is because we are using a simple model, and have put non-informative priors on our parameters.

One way to visualize the variability in our estimation of the regression line is to plot multiple estimates from the posterior.

library(gdata)

plot(y ~ x, pch = 20)

for (i in 1:500) {
 abline(posterior$alpha[i], posterior$beta[i], col = "gray", lty = 1)
}

abline(mean(posterior$alpha), mean(posterior$beta), col = 6, lw = 2)

####################################################################

Let’s try again, but now with more informative priors for the relationship between sea ice and time. We’re going to use normal priors with small standard deviations. If we were to use normal priors with very large standard deviations (say 1000, or 10,000), they would act very similarly to uniform priors.

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 alpha ~ normal(10, 0.1);
 beta ~ normal(1, 0.1);
 y ~ normal(alpha + x * beta , sigma);
}

generated quantities {}",

"stan_model2.stan")

stan_model2 <- "stan_model2.stan"

fit2 <- stan(stan_model2, data = stan_data, warmup = 500, iter = 1000, chains = 4, cores = 2, thin = 1)

posterior2 <- extract(fit2)

plot(y ~ x, pch = 20)

#abline(alpha, beta, col = 4, lty = 2, lw = 2)
abline(mean(posterior2$alpha), mean(posterior2$beta), col = 3, lw = 2)
abline(mean(posterior$alpha), mean(posterior$beta), col = 36, lw = 3)

This is a common issue in Bayesian modelling, if your prior distributions are very narrow and yet don’t fit your understanding of the system or the distribution of your data, you could run models that do not meaningfully explain variation in your data. However, that isn’t to say that you shouldn’t choose somewhat informative priors, you do want to use previous analyses and understanding of your study system inform your model priors and design. You just need to think carefully about each modelling decision you make!

Convergence Diagnostics For traceplots, we can view them directly from the posterior:

plot(posterior$alpha, type = "l")

plot(posterior$beta, type = "l")

plot(posterior$sigma, type = "l")

traceplot(stan_model1.fit)

Example of poor convergence

Try running a model for only 50 iterations and check the traceplots.

fit_bad <- stan(stan_model1, data = stan_data, warmup = 25, iter = 50, chains = 4, cores = 2, thin = 1)
## Warning: There were 5 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## Warning: Examine the pairs() plot to diagnose sampling problems
posterior_bad <- extract(fit_bad)
plot(posterior_bad$alpha, type = "l")

plot(posterior_bad$beta, type = "l")

plot(posterior_bad$sigma, type = "l")

Parameter summaries

We can also get summaries of the parameters through the posterior directly. Let’s also plot the non-Bayesian linear model values to make sure our model is doing what we think it is…

par(mfrow = c(1,3))

plot(density(posterior$alpha), main = "Alpha")
abline(v = lm_alpha, col = 4, lty = 2)

plot(density(posterior$beta), main = "Beta")
abline(v = lm_beta, col = 4, lty = 2)

plot(density(posterior$sigma), main = "Sigma")
abline(v = lm_sigma, col = 4, lty = 2)

From the posterior we can directly calculate the probability of any parameter being over or under a certain value of interest.

#Probablility that beta is >0:

sum(posterior$beta>0)/length(posterior$beta)
## [1] 0
#Probablility that beta is >0.2:

sum(posterior$beta>0.2)/length(posterior$beta)
## [1] 0

Posterior Predictive Checks

For prediction and as another form of model diagnostic, Stan can use random number generators to generate predicted values for each data point, at each iteration. This way we can generate predictions that also represent the uncertainties in our model and our data generation process. We generate these using the Generated Quantities block. This block can be used to get any other information we want about the posterior, or make predictions for new data.

Note that vectorization is not supported in the GQ (generated quantities) block, so we have to put it in a loop. But since this is compiled to C++, loops are actually quite fast and Stan only evaluates the GQ block once per iteration, so it won’t add too much time to your sampling. Typically, the data generating functions will be the distributions you used in the model block but with an _rng suffix. (Double-check in the Stan manual to see which sampling statements have corresponding rng functions already coded up.)

write("// Stan model for simple linear regression

data {
 int < lower = 1 > N; // Sample size
 vector[N] x; // Predictor
 vector[N] y; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 y ~ normal(x * beta + alpha, sigma);
}

generated quantities {
 real y_rep[N];

 for (n in 1:N) {
 y_rep[n] = normal_rng(x[n] * beta + alpha, sigma);
 }

}",

"stan_model2_GQ.stan")

stan_model2_GQ <- "stan_model2_GQ.stan"

fit3 <- stan(stan_model2_GQ, data = stan_data, iter = 1000, chains = 4, cores = 2, thin = 1)

#Extracting the y_rep values from posterior.
#Each row is an iteration (single posterior estimate) from the model.
y_rep <- as.matrix(fit3, pars = "y_rep")
dim(y_rep)
## [1] 2000   39
#Comparing density of y with densities of y over 200 posterior draws.

bayesplot::ppc_dens_overlay(y, y_rep[1:200, ])

#We can also use this to compare estimates of summary statistics.
bayesplot::ppc_stat(y = y, yrep = y_rep, stat = "mean")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#We can investigate mean posterior prediction per datapoint vs 
#the observed value for each datapoint (default line is 1:1)

bayesplot::ppc_scatter_avg(y = y, yrep = y_rep)

#Here is a list of currently available plots (bayesplot 1.2):
#available_ppc()
#can change the colour scheme in bayesplot too:
bayesplot::color_scheme_view(c("blue", "gray", "green", "pink", "purple",
 "red","teal","yellow"))

#And you can even mix them:

color_scheme_view("mix-blue-red")

#You can set color schemes with:

color_scheme_set("blue")

0.9.3 rats example

data {
  int<lower=0> N;
  int<lower=0> T;
  real x[T];
  real y[N,T];
  real xbar;
}
parameters {
  real alpha[N];
  real beta[N];

  real mu_alpha;
  real mu_beta;          // beta.c in original bugs model

  real<lower=0> sigmasq_y;
  real<lower=0> sigmasq_alpha;
  real<lower=0> sigmasq_beta;
}
transformed parameters {
  real<lower=0> sigma_y;       // sigma in original bugs model
  real<lower=0> sigma_alpha;
  real<lower=0> sigma_beta;

  sigma_y = sqrt(sigmasq_y);
  sigma_alpha = sqrt(sigmasq_alpha);
  sigma_beta =  sqrt(sigmasq_beta);
}
model {
  mu_alpha ~ normal(0, 100);
  mu_beta ~ normal(0, 100);
  sigmasq_y ~ inv_gamma(0.001, 0.001);
  sigmasq_alpha ~ inv_gamma(0.001, 0.001);
  sigmasq_beta ~ inv_gamma(0.001, 0.001);
  alpha ~ normal(mu_alpha, sigma_alpha); // vectorized
  beta ~ normal(mu_beta, sigma_beta);  // vectorized
  for (n in 1:N)
    for (t in 1:T) 
      y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), sigma_y);

}
generated quantities {
  real alpha0;
  alpha0 = mu_alpha - xbar * mu_beta;
}
library(tidyverse)
df <- read_delim("rats.txt", delim = " ")
## Parsed with column specification:
## cols(
##   day8 = col_double(),
##   day15 = col_double(),
##   day22 = col_double(),
##   day29 = col_double(),
##   day36 = col_double()
## )
df
## # A tibble: 30 x 5
##     day8 day15 day22 day29 day36
##    <dbl> <dbl> <dbl> <dbl> <dbl>
##  1   151   199   246   283   320
##  2   145   199   249   293   354
##  3   147   214   263   312   328
##  4   155   200   237   272   297
##  5   135   188   230   280   323
##  6   159   210   252   298   331
##  7   141   189   231   275   305
##  8   159   201   248   297   338
##  9   177   236   285   350   376
## 10   134   182   220   260   296
## # … with 20 more rows
y <- as.matrix(df)
y
##       day8 day15 day22 day29 day36
##  [1,]  151   199   246   283   320
##  [2,]  145   199   249   293   354
##  [3,]  147   214   263   312   328
##  [4,]  155   200   237   272   297
##  [5,]  135   188   230   280   323
##  [6,]  159   210   252   298   331
##  [7,]  141   189   231   275   305
##  [8,]  159   201   248   297   338
##  [9,]  177   236   285   350   376
## [10,]  134   182   220   260   296
## [11,]  160   208   261   313   352
## [12,]  143   188   220   273   314
## [13,]  154   200   244   289   325
## [14,]  171   221   270   326   358
## [15,]  163   216   242   281   312
## [16,]  160   207   248   288   324
## [17,]  142   187   234   280   316
## [18,]  156   203   243   283   317
## [19,]  157   212   259   307   336
## [20,]  152   203   246   286   321
## [21,]  154   205   253   298   334
## [22,]  139   190   225   267   302
## [23,]  146   191   229   272   302
## [24,]  157   211   250   285   323
## [25,]  132   185   237   286   331
## [26,]  160   207   257   303   345
## [27,]  169   216   261   295   333
## [28,]  157   205   248   289   316
## [29,]  137   180   219   258   291
## [30,]  153   200   244   286   324
x <- c(8,15,22,29,36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
library(rstan)
options(mc.cores = parallel::detectCores())
#sample from that posterior distribution 
rats_fit <- rstan::sampling(rats, 
                     data = list(y,x,xbar,N,T)) 
rstan::summary(rats_fit)
## $summary
##                       mean     se_mean         sd         2.5%
## alpha[1]       239.9101176 0.032107296  2.6578237  234.5488751
## alpha[2]       247.7710919 0.029412775  2.6246964  242.6236101
## alpha[3]       252.4522863 0.031460374  2.6109330  247.2652729
## alpha[4]       232.5657126 0.030989675  2.7144837  227.1594045
## alpha[5]       231.6309075 0.034440381  2.7045632  226.3269739
## alpha[6]       249.7481589 0.030889327  2.6572009  244.4669857
## alpha[7]       228.7758177 0.032963901  2.7019954  223.4477749
## alpha[8]       248.3437472 0.033694844  2.7709166  243.0794437
## alpha[9]       283.3338365 0.034574598  2.7318716  277.8795001
## alpha[10]      219.2564254 0.036961968  2.6462226  214.0145444
## alpha[11]      258.2312268 0.036763933  2.7278808  252.8447814
## alpha[12]      228.1444671 0.037685651  2.7347577  222.7936056
## alpha[13]      242.3802764 0.034283086  2.7192027  237.1123343
## alpha[14]      268.2839619 0.035878255  2.7584002  262.7762843
## alpha[15]      242.8368719 0.033333623  2.6940760  237.7157202
## alpha[16]      245.3301513 0.031967736  2.7005467  239.9721467
## alpha[17]      232.1730819 0.036344901  2.7550520  226.7701743
## alpha[18]      240.5326904 0.035837571  2.7692615  235.0787279
## alpha[19]      253.7550682 0.034482225  2.6786253  248.4006413
## alpha[20]      241.6498238 0.034534730  2.7206198  236.2303313
## alpha[21]      248.6097520 0.033301156  2.7067765  243.4666249
## alpha[22]      225.2224625 0.036639060  2.6407778  220.1035482
## alpha[23]      228.5369983 0.037710436  2.6966543  223.1865699
## alpha[24]      245.1108964 0.032851918  2.6535119  239.9802538
## alpha[25]      234.5014763 0.030387887  2.7142259  229.0450572
## alpha[26]      253.9817892 0.034379196  2.7662874  248.4171417
## alpha[27]      254.3184138 0.033148950  2.6496400  249.2526166
## alpha[28]      242.9907623 0.034766967  2.6758857  237.5824858
## alpha[29]      217.9332278 0.035070957  2.7041637  212.7331971
## alpha[30]      241.4534414 0.034657623  2.7229125  236.2199885
## beta[1]          6.0644501 0.003317410  0.2450428    5.5844573
## beta[2]          7.0528794 0.003605541  0.2557732    6.5299716
## beta[3]          6.4840340 0.003160061  0.2529363    5.9827008
## beta[4]          5.3421506 0.003280252  0.2528233    4.8292343
## beta[5]          6.5693428 0.003064760  0.2402844    6.0875833
## beta[6]          6.1741715 0.003276169  0.2432543    5.6955275
## beta[7]          5.9754878 0.003019117  0.2421691    5.4946117
## beta[8]          6.4204992 0.003119430  0.2468635    5.9338214
## beta[9]          7.0527735 0.003306984  0.2601278    6.5253891
## beta[10]         5.8549368 0.002974570  0.2434404    5.3741981
## beta[11]         6.7986044 0.003260923  0.2556692    6.2833233
## beta[12]         6.1200303 0.002905360  0.2418928    5.6448467
## beta[13]         6.1684915 0.003273787  0.2495856    5.6505664
## beta[14]         6.6907145 0.002901221  0.2580016    6.1858630
## beta[15]         5.4163624 0.003153302  0.2442228    4.9372132
## beta[16]         5.9217063 0.003038135  0.2404501    5.4428470
## beta[17]         6.2763483 0.003103078  0.2465650    5.7935908
## beta[18]         5.8473081 0.002921502  0.2426703    5.3645122
## beta[19]         6.4037304 0.003107552  0.2449235    5.9266222
## beta[20]         6.0547395 0.003046638  0.2448963    5.5626759
## beta[21]         6.4024657 0.003161328  0.2449800    5.9147997
## beta[22]         5.8627274 0.003171593  0.2435659    5.3924864
## beta[23]         5.7456716 0.003123035  0.2516111    5.2683213
## beta[24]         5.8917624 0.002771738  0.2472508    5.3973264
## beta[25]         6.9061627 0.003119167  0.2457827    6.4230444
## beta[26]         6.5481125 0.003068555  0.2427410    6.0677120
## beta[27]         5.9031704 0.002825252  0.2389190    5.4280174
## beta[28]         5.8487791 0.002891054  0.2440778    5.3605305
## beta[29]         5.6742810 0.003277911  0.2464774    5.1871659
## beta[30]         6.1342841 0.002726781  0.2374156    5.6732897
## mu_alpha       242.4769641 0.036647467  2.7496143  236.9904708
## mu_beta          6.1873138 0.001650494  0.1083687    5.9689339
## sigmasq_y       37.4509648 0.111887922  5.8444217   27.3319211
## sigmasq_alpha  217.8323501 1.135202289 64.6337673  124.8698817
## sigmasq_beta     0.2735714 0.001689988  0.1000803    0.1315568
## sigma_y          6.1015115 0.009067297  0.4717821    5.2279940
## sigma_alpha     14.6086201 0.034923241  2.1027775   11.1745193
## sigma_beta       0.5149685 0.001535668  0.0915475    0.3627077
## alpha0         106.3560610 0.052391111  3.6219112   99.0589096
## lp__          -438.6669450 0.216886901  7.0499385 -453.4130373
##                        25%          50%          75%        97.5%    n_eff
## alpha[1]       238.1675377  239.9333958  241.6984956  245.0355233 6852.434
## alpha[2]       246.0733276  247.7323203  249.5515871  253.0381621 7963.173
## alpha[3]       250.7041369  252.4513488  254.2281914  257.5878145 6887.533
## alpha[4]       230.7353875  232.5541032  234.3574918  237.8531382 7672.562
## alpha[5]       229.7891779  231.5972158  233.4234115  236.8976075 6166.779
## alpha[6]       247.9537566  249.7786087  251.5766478  254.8018159 7400.003
## alpha[7]       226.9337347  228.7873045  230.6294379  234.0227400 6718.804
## alpha[8]       246.4098075  248.2996921  250.2499589  253.9513876 6762.698
## alpha[9]       281.5543810  283.3615450  285.2205777  288.5110862 6243.186
## alpha[10]      217.4248965  219.2374032  221.0676174  224.4912249 5125.575
## alpha[11]      256.4146007  258.2244878  260.0467161  263.6088276 5505.628
## alpha[12]      226.3276403  228.2052378  229.9760624  233.5149403 5266.058
## alpha[13]      240.5610480  242.3535576  244.2290557  247.6511055 6291.053
## alpha[14]      266.4408861  268.3543095  270.1169785  273.6909386 5910.877
## alpha[15]      241.0688565  242.8195377  244.6080161  248.1623050 6532.127
## alpha[16]      243.5343836  245.3560795  247.1501472  250.4967635 7136.407
## alpha[17]      230.2873536  232.2097496  234.0473887  237.4544068 5746.093
## alpha[18]      238.7007733  240.5132073  242.3822496  246.0491638 5971.051
## alpha[19]      251.9807564  253.7744143  255.5407144  259.0320223 6034.390
## alpha[20]      239.8380285  241.6809851  243.4791093  247.0318862 6206.169
## alpha[21]      246.7655637  248.5731640  250.4145555  254.1298452 6606.724
## alpha[22]      223.4544631  225.2162203  226.9496564  230.3653646 5194.875
## alpha[23]      226.7496189  228.5775870  230.3223556  233.7823100 5113.608
## alpha[24]      243.3497722  245.1177653  246.8428643  250.3437867 6524.100
## alpha[25]      232.7671886  234.5388923  236.2405042  239.7228800 7977.944
## alpha[26]      252.1372000  253.9935830  255.8300567  259.4362636 6474.454
## alpha[27]      252.4940461  254.3519698  256.1306437  259.5492879 6389.019
## alpha[28]      241.1898075  242.9504615  244.7909839  248.3932803 5923.815
## alpha[29]      216.1015136  217.8699391  219.7309097  223.3787407 5945.259
## alpha[30]      239.6092495  241.4320058  243.2925863  246.7302783 6172.624
## beta[1]          5.9016868    6.0602342    6.2295013    6.5420962 5456.140
## beta[2]          6.8826410    7.0528893    7.2319203    7.5503089 5032.332
## beta[3]          6.3180118    6.4857643    6.6544347    6.9711355 6406.654
## beta[4]          5.1757242    5.3423518    5.5111674    5.8414451 5940.455
## beta[5]          6.4076781    6.5706441    6.7288095    7.0411539 6146.929
## beta[6]          6.0091246    6.1776375    6.3427587    6.6495631 5513.005
## beta[7]          5.8158684    5.9753670    6.1360390    6.4496104 6433.946
## beta[8]          6.2573140    6.4210072    6.5842162    6.9064715 6262.721
## beta[9]          6.8829068    7.0492251    7.2239928    7.5697179 6187.417
## beta[10]         5.6963592    5.8553535    6.0152041    6.3303438 6697.871
## beta[11]         6.6320326    6.8038259    6.9667485    7.3088905 6147.175
## beta[12]         5.9624477    6.1203819    6.2827176    6.5943635 6931.800
## beta[13]         6.0050765    6.1685567    6.3315431    6.6565623 5812.169
## beta[14]         6.5156479    6.6855874    6.8674784    7.1907583 7908.303
## beta[15]         5.2558817    5.4176733    5.5789920    5.9111793 5998.478
## beta[16]         5.7599369    5.9215235    6.0876286    6.3764558 6263.768
## beta[17]         6.1162209    6.2763330    6.4426962    6.7593673 6313.605
## beta[18]         5.6871604    5.8473301    6.0057197    6.3329377 6899.549
## beta[19]         6.2328684    6.4047654    6.5708075    6.8761905 6211.893
## beta[20]         5.8913325    6.0545078    6.2197273    6.5378323 6461.342
## beta[21]         6.2385311    6.4065739    6.5686046    6.8716331 6005.127
## beta[22]         5.6945705    5.8644918    6.0289380    6.3420863 5897.635
## beta[23]         5.5691256    5.7417880    5.9176685    6.2343153 6490.916
## beta[24]         5.7303897    5.8956260    6.0487064    6.3805677 7957.397
## beta[25]         6.7425226    6.9051803    7.0704856    7.3998105 6209.051
## beta[26]         6.3902110    6.5469365    6.7130224    7.0243083 6257.755
## beta[27]         5.7390866    5.9034882    6.0647555    6.3668808 7151.333
## beta[28]         5.6875684    5.8456505    6.0137796    6.3366710 7127.613
## beta[29]         5.5104588    5.6764913    5.8348260    6.1467949 5654.051
## beta[30]         5.9759472    6.1376583    6.2936226    6.6079329 7580.853
## mu_alpha       240.6768502  242.5065350  244.3402045  247.6967971 5629.316
## mu_beta          6.1156840    6.1904349    6.2572214    6.4017220 4311.020
## sigmasq_y       33.2734113   37.0175766   41.0758986   50.2616686 2728.454
## sigmasq_alpha  171.5328956  206.8671150  251.6332633  376.1496098 3241.697
## sigmasq_beta     0.2021241    0.2569905    0.3234950    0.5170319 3506.953
## sigma_y          5.7683110    6.0842071    6.4090482    7.0895462 2707.243
## sigma_alpha     13.0970568   14.3828757   15.8629525   19.3945767 3625.414
## sigma_beta       0.4495821    0.5069423    0.5687662    0.7190493 3553.844
## alpha0         103.9324694  106.3789046  108.7889424  113.3425357 4779.257
## lp__          -443.1543728 -438.4595874 -433.7743000 -425.7954482 1056.584
##                    Rhat
## alpha[1]      1.0000846
## alpha[2]      0.9992497
## alpha[3]      0.9996211
## alpha[4]      0.9996334
## alpha[5]      0.9992586
## alpha[6]      0.9996418
## alpha[7]      0.9993552
## alpha[8]      1.0001294
## alpha[9]      0.9994213
## alpha[10]     1.0003155
## alpha[11]     0.9996596
## alpha[12]     0.9995462
## alpha[13]     0.9992462
## alpha[14]     0.9995434
## alpha[15]     0.9996318
## alpha[16]     0.9991470
## alpha[17]     0.9993781
## alpha[18]     0.9995428
## alpha[19]     0.9991064
## alpha[20]     0.9998488
## alpha[21]     1.0001483
## alpha[22]     1.0000499
## alpha[23]     0.9999544
## alpha[24]     0.9993894
## alpha[25]     1.0000089
## alpha[26]     0.9994473
## alpha[27]     0.9999175
## alpha[28]     0.9991887
## alpha[29]     0.9996949
## alpha[30]     0.9993756
## beta[1]       0.9997811
## beta[2]       0.9993809
## beta[3]       0.9992219
## beta[4]       0.9995006
## beta[5]       0.9995295
## beta[6]       1.0000908
## beta[7]       0.9995685
## beta[8]       0.9994455
## beta[9]       0.9992980
## beta[10]      0.9995895
## beta[11]      0.9994011
## beta[12]      0.9994513
## beta[13]      1.0001090
## beta[14]      0.9997021
## beta[15]      0.9995223
## beta[16]      0.9993061
## beta[17]      0.9994898
## beta[18]      0.9993602
## beta[19]      1.0001532
## beta[20]      0.9997914
## beta[21]      0.9995996
## beta[22]      0.9999703
## beta[23]      0.9992240
## beta[24]      0.9991987
## beta[25]      0.9996382
## beta[26]      0.9995108
## beta[27]      0.9998728
## beta[28]      0.9998624
## beta[29]      0.9993781
## beta[30]      0.9992577
## mu_alpha      0.9996557
## mu_beta       0.9994931
## sigmasq_y     1.0005936
## sigmasq_alpha 1.0015092
## sigmasq_beta  1.0001771
## sigma_y       1.0005634
## sigma_alpha   1.0012839
## sigma_beta    1.0001317
## alpha0        0.9996848
## lp__          1.0021535
## 
## $c_summary
## , , chains = chain:1
## 
##                stats
## parameter               mean          sd         2.5%          25%
##   alpha[1]       239.7972483  2.60639912  234.7125902  238.0280946
##   alpha[2]       247.8251074  2.64178150  242.6955335  246.0865869
##   alpha[3]       252.4075042  2.42149759  247.8130702  250.6569773
##   alpha[4]       232.5820993  2.78701190  227.1576092  230.7455688
##   alpha[5]       231.6376926  2.84620844  226.1038192  229.6695463
##   alpha[6]       249.8192304  2.53488114  245.0315400  248.0620879
##   alpha[7]       228.7994238  2.63838802  223.4539937  227.0364576
##   alpha[8]       248.5014452  2.74863684  243.1739216  246.6413974
##   alpha[9]       283.2586783  2.81259040  277.8128621  281.4368390
##   alpha[10]      219.1674360  2.67444720  214.0436032  217.2310065
##   alpha[11]      258.1701723  2.64754811  253.0265912  256.3253384
##   alpha[12]      228.1554644  2.79487764  222.7650150  226.3199139
##   alpha[13]      242.3428627  2.65700499  237.0236869  240.5946878
##   alpha[14]      268.2743853  2.79573460  262.7344352  266.4646927
##   alpha[15]      242.9084497  2.59779595  237.8332553  241.2410825
##   alpha[16]      245.3464185  2.77616064  239.8875378  243.4896350
##   alpha[17]      232.1369169  2.67040252  226.7696592  230.2940502
##   alpha[18]      240.4450581  2.83306438  234.8916814  238.4722564
##   alpha[19]      253.7630544  2.67810152  248.4073626  251.9712172
##   alpha[20]      241.5329743  2.63624174  236.5329043  239.6132601
##   alpha[21]      248.7009847  2.79779759  243.2702018  246.8070504
##   alpha[22]      225.1551859  2.55138412  220.0214701  223.5254486
##   alpha[23]      228.4582150  2.66217730  223.5799146  226.6417695
##   alpha[24]      245.0845830  2.64331907  240.1111242  243.3219981
##   alpha[25]      234.4808369  2.65146943  228.9687418  232.8232873
##   alpha[26]      253.9711000  2.71312087  248.4631016  252.2982053
##   alpha[27]      254.4157853  2.51900937  249.3283539  252.6990020
##   alpha[28]      243.0326897  2.50076114  238.0271618  241.3584012
##   alpha[29]      217.9333756  2.63575985  212.8333935  216.1945747
##   alpha[30]      241.4935969  2.67247409  236.3254160  239.6412567
##   beta[1]          6.0725365  0.25375758    5.5751724    5.9005951
##   beta[2]          7.0572707  0.26576382    6.5272554    6.8768346
##   beta[3]          6.4785965  0.25848103    5.9545301    6.3103362
##   beta[4]          5.3378774  0.25591881    4.8546500    5.1632625
##   beta[5]          6.5660464  0.23790468    6.1005313    6.4013382
##   beta[6]          6.1681612  0.25476636    5.6588987    5.9948806
##   beta[7]          5.9794942  0.23094359    5.5484815    5.8200887
##   beta[8]          6.4250619  0.23607641    5.9714962    6.2638567
##   beta[9]          7.0479970  0.26126136    6.5210483    6.8780313
##   beta[10]         5.8554059  0.23998411    5.3694576    5.6999254
##   beta[11]         6.8014372  0.25903763    6.2636761    6.6329094
##   beta[12]         6.1278073  0.24319484    5.6680165    5.9540386
##   beta[13]         6.1621114  0.23798081    5.6718732    6.0039676
##   beta[14]         6.6924534  0.25517657    6.1751998    6.5243266
##   beta[15]         5.4117894  0.25487922    4.9118082    5.2518834
##   beta[16]         5.9153507  0.24062944    5.4173267    5.7554914
##   beta[17]         6.2856552  0.23759743    5.8219599    6.1238942
##   beta[18]         5.8436079  0.23358626    5.3756507    5.6839407
##   beta[19]         6.3917618  0.24054249    5.9459660    6.2242411
##   beta[20]         6.0612853  0.24378011    5.5561171    5.8987808
##   beta[21]         6.4050348  0.23185595    5.9690453    6.2434648
##   beta[22]         5.8692426  0.24153703    5.3892752    5.7122298
##   beta[23]         5.7450932  0.25431074    5.2639867    5.5627092
##   beta[24]         5.8953821  0.24383241    5.4019218    5.7432739
##   beta[25]         6.9151386  0.24303343    6.4107486    6.7615367
##   beta[26]         6.5509966  0.24114268    6.0755573    6.3911304
##   beta[27]         5.9091451  0.23952422    5.4249514    5.7354093
##   beta[28]         5.8409380  0.24879376    5.3467413    5.6777638
##   beta[29]         5.6679922  0.25588099    5.1568610    5.4990600
##   beta[30]         6.1339759  0.22812373    5.6922933    5.9931439
##   mu_alpha       242.4409566  2.72575406  237.0569237  240.5731728
##   mu_beta          6.1888002  0.10869891    5.9644545    6.1195778
##   sigmasq_y       37.2086601  5.75626152   27.9376150   33.0239414
##   sigmasq_alpha  220.0866857 65.32654435  126.9929368  172.8760351
##   sigmasq_beta     0.2779605  0.10214142    0.1356473    0.2026671
##   sigma_y          6.0820942  0.46584020    5.2856021    5.7466461
##   sigma_alpha     14.6857275  2.10250193   11.2691142   13.1482331
##   sigma_beta       0.5190288  0.09261822    0.3683032    0.4501857
##   alpha0         106.2873529  3.59541168   99.1870589  103.8453352
##   lp__          -438.2808538  6.93621723 -452.9071952 -442.7341882
##                stats
## parameter                50%          75%        97.5%
##   alpha[1]       239.8667947  241.5703332  244.7670055
##   alpha[2]       247.7289901  249.5734675  253.1452331
##   alpha[3]       252.3328409  254.0046945  257.0674777
##   alpha[4]       232.5285156  234.3770633  238.0928150
##   alpha[5]       231.7535331  233.5433699  237.1264280
##   alpha[6]       249.8112129  251.5716195  254.5118904
##   alpha[7]       228.8476008  230.6032508  233.8660793
##   alpha[8]       248.3551000  250.3753860  254.0434512
##   alpha[9]       283.2527414  285.1978995  288.6836865
##   alpha[10]      219.1490321  221.1502726  224.2661844
##   alpha[11]      258.1900953  259.9643330  263.4089542
##   alpha[12]      228.1879610  229.9203585  233.6582897
##   alpha[13]      242.2845373  244.1705987  247.5685440
##   alpha[14]      268.3511802  270.0927333  273.8345933
##   alpha[15]      242.9757276  244.4742942  248.2689828
##   alpha[16]      245.4705045  247.2927467  250.5375328
##   alpha[17]      232.1878184  234.0007558  237.1797169
##   alpha[18]      240.4619626  242.4841965  245.6027199
##   alpha[19]      253.7804844  255.6277787  259.0771619
##   alpha[20]      241.5009428  243.3476970  246.8651783
##   alpha[21]      248.6055741  250.5588138  254.3670005
##   alpha[22]      225.1700510  226.8390154  229.8768042
##   alpha[23]      228.3636389  230.2150977  233.7860752
##   alpha[24]      245.0909185  246.8447215  250.0440502
##   alpha[25]      234.5708732  236.1869283  239.5833601
##   alpha[26]      253.9671944  255.7759811  259.4650977
##   alpha[27]      254.4697205  256.1536545  259.1550804
##   alpha[28]      243.0038281  244.7258425  248.0647944
##   alpha[29]      217.8881224  219.5746914  223.2441091
##   alpha[30]      241.4579842  243.2448924  246.5852838
##   beta[1]          6.0752151    6.2493067    6.5823282
##   beta[2]          7.0645709    7.2425927    7.5683913
##   beta[3]          6.4803484    6.6487742    6.9810153
##   beta[4]          5.3360573    5.5092392    5.8555246
##   beta[5]          6.5677169    6.7122138    7.0615434
##   beta[6]          6.1825373    6.3485085    6.6462179
##   beta[7]          5.9774969    6.1335177    6.4464679
##   beta[8]          6.4252280    6.5826515    6.9044311
##   beta[9]          7.0467226    7.2209842    7.5427927
##   beta[10]         5.8590136    6.0134038    6.3056013
##   beta[11]         6.8035609    6.9767230    7.3106500
##   beta[12]         6.1312030    6.2872747    6.5947122
##   beta[13]         6.1648139    6.3202458    6.6258229
##   beta[14]         6.6889807    6.8719880    7.1793983
##   beta[15]         5.4098748    5.5856729    5.9164090
##   beta[16]         5.9149716    6.0831396    6.3817143
##   beta[17]         6.2803634    6.4473034    6.7484575
##   beta[18]         5.8530618    5.9988373    6.2735124
##   beta[19]         6.3896452    6.5594383    6.8673249
##   beta[20]         6.0612798    6.2280326    6.5398203
##   beta[21]         6.4044263    6.5531506    6.8705243
##   beta[22]         5.8692256    6.0261199    6.3153091
##   beta[23]         5.7492518    5.9173904    6.2195379
##   beta[24]         5.8967558    6.0503862    6.3795245
##   beta[25]         6.9226011    7.0701602    7.3823138
##   beta[26]         6.5478469    6.7157076    7.0073336
##   beta[27]         5.9081767    6.0759034    6.3783507
##   beta[28]         5.8415502    6.0020651    6.3185916
##   beta[29]         5.6672162    5.8308601    6.1617017
##   beta[30]         6.1387495    6.2678595    6.5970509
##   mu_alpha       242.4736702  244.3411848  247.6314297
##   mu_beta          6.1911350    6.2561357    6.4063237
##   sigmasq_y       36.7446865   40.7849134   49.9108488
##   sigmasq_alpha  210.0478600  253.4239247  385.0172292
##   sigmasq_beta     0.2639536    0.3309054    0.5193950
##   sigma_y          6.0617396    6.3863067    7.0647606
##   sigma_alpha     14.4930278   15.9192941   19.6218359
##   sigma_beta       0.5137641    0.5752438    0.7206901
##   alpha0         106.2827742  108.7140118  112.9597328
##   lp__          -438.1609888 -433.2515897 -426.1289985
## 
## , , chains = chain:2
## 
##                stats
## parameter               mean          sd         2.5%          25%
##   alpha[1]       240.0747142  2.67480931  234.6014977  238.3510757
##   alpha[2]       247.7791483  2.48076409  242.8598396  246.1585024
##   alpha[3]       252.5400268  2.73344171  246.9316964  250.8053271
##   alpha[4]       232.5043000  2.73911279  227.0104892  230.6651318
##   alpha[5]       231.5724667  2.79956394  226.1547100  229.7264167
##   alpha[6]       249.6979044  2.66348224  244.2892738  247.9633249
##   alpha[7]       228.7417049  2.78991710  223.4282040  226.8823485
##   alpha[8]       248.3748856  2.82742957  242.7569469  246.2869964
##   alpha[9]       283.4186550  2.67413497  278.2628846  281.6506620
##   alpha[10]      219.4607838  2.61574465  214.5820323  217.7577664
##   alpha[11]      258.2763834  2.71785181  253.0133848  256.4410328
##   alpha[12]      228.0904648  2.70402805  222.7266779  226.3245153
##   alpha[13]      242.4535014  2.70159458  237.2739232  240.6666538
##   alpha[14]      268.3856271  2.84523616  262.8407007  266.4889397
##   alpha[15]      242.7291418  2.66002406  237.7145691  241.0109200
##   alpha[16]      245.3076640  2.70145230  239.9451239  243.5293278
##   alpha[17]      232.1204437  2.65399778  227.1262676  230.2473892
##   alpha[18]      240.5630063  2.71429868  235.2813937  238.6986005
##   alpha[19]      253.7622650  2.80493841  248.1296828  251.7788658
##   alpha[20]      241.6292937  2.74220110  236.1444055  239.8148044
##   alpha[21]      248.5653116  2.65089624  243.5092800  246.8583196
##   alpha[22]      225.2495660  2.62349610  220.4590882  223.4004887
##   alpha[23]      228.5922810  2.83136168  223.1915052  226.6264426
##   alpha[24]      245.1382897  2.59282809  240.1975897  243.4616159
##   alpha[25]      234.6212015  2.72561285  229.3240090  232.8437234
##   alpha[26]      254.0276301  2.85061782  248.7552747  251.9345470
##   alpha[27]      254.2933415  2.73820219  249.2976326  252.1929444
##   alpha[28]      242.9818848  2.72103431  237.7555241  241.0687728
##   alpha[29]      217.9564911  2.77248359  212.8217920  215.9531638
##   alpha[30]      241.4062730  2.78460016  235.9520482  239.6014886
##   beta[1]          6.0570757  0.23995987    5.5592126    5.9031738
##   beta[2]          7.0455736  0.24823125    6.5486288    6.8808357
##   beta[3]          6.4872947  0.24987883    5.9958932    6.3317133
##   beta[4]          5.3374653  0.25112179    4.8184346    5.1709624
##   beta[5]          6.5673969  0.23806505    6.0984647    6.4075377
##   beta[6]          6.1782835  0.24368123    5.6953206    6.0084617
##   beta[7]          5.9818976  0.24710363    5.4850908    5.8212416
##   beta[8]          6.4203687  0.25324816    5.9238352    6.2552395
##   beta[9]          7.0604632  0.25859440    6.5432910    6.8856921
##   beta[10]         5.8493531  0.25098768    5.3367616    5.6847870
##   beta[11]         6.7959306  0.24151337    6.3082028    6.6425787
##   beta[12]         6.1163597  0.23988947    5.6477409    5.9607595
##   beta[13]         6.1630874  0.24623088    5.6526444    6.0065985
##   beta[14]         6.6914393  0.26509687    6.1824109    6.4991512
##   beta[15]         5.4167923  0.23699726    4.9668301    5.2661803
##   beta[16]         5.9246342  0.23004250    5.4726913    5.7701535
##   beta[17]         6.2739346  0.24434763    5.7928751    6.1165333
##   beta[18]         5.8491454  0.25495552    5.3555984    5.6817647
##   beta[19]         6.4124240  0.24349639    5.9645970    6.2416357
##   beta[20]         6.0550317  0.24775494    5.5418660    5.8954408
##   beta[21]         6.3931212  0.24361230    5.9003630    6.2406859
##   beta[22]         5.8615080  0.25233071    5.3859543    5.6697783
##   beta[23]         5.7491811  0.24152150    5.2846098    5.5801701
##   beta[24]         5.8915084  0.23950650    5.4132766    5.7417108
##   beta[25]         6.9023390  0.24699475    6.4208513    6.7349543
##   beta[26]         6.5462247  0.23979237    6.0720912    6.3912508
##   beta[27]         5.8975327  0.24375454    5.4440496    5.7367358
##   beta[28]         5.8544269  0.24488914    5.3743382    5.6923970
##   beta[29]         5.6730817  0.23902089    5.2002456    5.5119920
##   beta[30]         6.1310352  0.24644089    5.6468563    5.9620273
##   mu_alpha       242.4516883  2.91564888  236.5638855  240.5695985
##   mu_beta          6.1864800  0.10850239    5.9719642    6.1142867
##   sigmasq_y       37.5397658  5.77165671   26.7517692   33.5337760
##   sigmasq_alpha  221.7872129 70.26058358  123.6443747  171.6186410
##   sigmasq_beta     0.2754960  0.09896423    0.1314987    0.2026667
##   sigma_y          6.1090332  0.46872082    5.1722112    5.7908354
##   sigma_alpha     14.7202996  2.25944628   11.1195484   13.1003297
##   sigma_beta       0.5169419  0.09096892    0.3626275    0.4501852
##   alpha0         106.3491284  3.69921527   98.8126777  103.9975807
##   lp__          -438.8659929  7.16304424 -455.2337088 -443.4893529
##                stats
## parameter                50%          75%        97.5%
##   alpha[1]       240.1610520  241.7474582  245.3130745
##   alpha[2]       247.7859673  249.4398650  252.9181520
##   alpha[3]       252.6512448  254.4022858  257.8744096
##   alpha[4]       232.4895592  234.2154667  238.0307830
##   alpha[5]       231.5143614  233.4810615  237.3391664
##   alpha[6]       249.7116914  251.4782609  254.6445344
##   alpha[7]       228.7621856  230.6846605  234.2117680
##   alpha[8]       248.4956397  250.4086693  253.9158397
##   alpha[9]       283.3664683  285.3723334  288.5136279
##   alpha[10]      219.3347147  221.2365367  224.7189544
##   alpha[11]      258.2335374  260.1207969  263.4907740
##   alpha[12]      228.2126934  229.9389770  233.2161719
##   alpha[13]      242.4077600  244.2092457  247.7374655
##   alpha[14]      268.3807680  270.2517317  273.8917225
##   alpha[15]      242.6550634  244.5996129  247.8651141
##   alpha[16]      245.3410139  247.0739714  250.3242025
##   alpha[17]      232.0843693  233.9659784  237.2120838
##   alpha[18]      240.5114448  242.4961966  246.0306898
##   alpha[19]      253.8131664  255.5839424  259.2131635
##   alpha[20]      241.6368721  243.5353991  247.2077761
##   alpha[21]      248.5596940  250.3542836  253.4632007
##   alpha[22]      225.1339228  227.0613126  230.3345192
##   alpha[23]      228.7652007  230.5233440  233.7702358
##   alpha[24]      245.1432921  246.7210075  250.3509015
##   alpha[25]      234.6890684  236.4008436  239.7644741
##   alpha[26]      253.9966584  256.0318602  259.1694514
##   alpha[27]      254.2857658  256.2760987  259.6261167
##   alpha[28]      242.9667833  244.8326251  248.6275956
##   alpha[29]      217.8886642  219.9075268  223.7470454
##   alpha[30]      241.3483505  243.3118613  246.5809743
##   beta[1]          6.0574730    6.2188376    6.5209861
##   beta[2]          7.0366741    7.2100928    7.5447620
##   beta[3]          6.4893272    6.6537480    6.9565601
##   beta[4]          5.3413451    5.5115431    5.8213470
##   beta[5]          6.5620287    6.7306732    7.0246114
##   beta[6]          6.1846659    6.3488964    6.6354444
##   beta[7]          5.9792160    6.1429901    6.4750550
##   beta[8]          6.4290286    6.5825511    6.8986123
##   beta[9]          7.0561164    7.2312302    7.5773047
##   beta[10]         5.8499079    6.0060804    6.3662178
##   beta[11]         6.8022274    6.9474881    7.2997560
##   beta[12]         6.1164113    6.2832467    6.5691386
##   beta[13]         6.1658161    6.3144293    6.6566333
##   beta[14]         6.6930047    6.8726290    7.1978639
##   beta[15]         5.4104574    5.5675423    5.9054396
##   beta[16]         5.9223063    6.0792712    6.3737321
##   beta[17]         6.2783930    6.4368432    6.7621670
##   beta[18]         5.8465627    6.0183761    6.3533501
##   beta[19]         6.4079061    6.5856453    6.8731707
##   beta[20]         6.0619902    6.2106457    6.5224253
##   beta[21]         6.3977751    6.5564404    6.8495474
##   beta[22]         5.8563616    6.0432677    6.3627721
##   beta[23]         5.7498557    5.9144709    6.2191615
##   beta[24]         5.8958496    6.0408314    6.3692476
##   beta[25]         6.8971236    7.0724173    7.3810687
##   beta[26]         6.5425162    6.7025793    7.0243543
##   beta[27]         5.8968628    6.0601252    6.3730498
##   beta[28]         5.8512035    6.0193571    6.3423237
##   beta[29]         5.6722369    5.8355217    6.1295456
##   beta[30]         6.1397007    6.2937324    6.6235798
##   mu_alpha       242.4641938  244.4273035  247.8313812
##   mu_beta          6.1887881    6.2585127    6.3987412
##   sigmasq_y       37.1891048   41.1386489   50.0381626
##   sigmasq_alpha  207.7786588  256.6538142  398.3097495
##   sigmasq_beta     0.2586505    0.3282218    0.4925724
##   sigma_y          6.0982870    6.4139417    7.0737654
##   sigma_alpha     14.4145293   16.0204179   19.9576989
##   sigma_beta       0.5085770    0.5729063    0.7018350
##   alpha0         106.3449729  108.8658653  113.2457998
##   lp__          -438.4862993 -433.9871375 -425.5198251
## 
## , , chains = chain:3
## 
##                stats
## parameter               mean          sd         2.5%          25%
##   alpha[1]       239.8817026  2.69856802  234.5003184  238.0937825
##   alpha[2]       247.6981732  2.74071929  242.1202433  245.9699330
##   alpha[3]       252.4725536  2.61396199  247.2347767  250.7270517
##   alpha[4]       232.6650683  2.58106078  227.3071115  230.9328126
##   alpha[5]       231.6676142  2.62784112  226.5553436  229.8802622
##   alpha[6]       249.6630372  2.65022329  244.5701236  247.9245898
##   alpha[7]       228.7975073  2.83303105  223.1622942  226.7577078
##   alpha[8]       248.2318881  2.77037261  242.9606118  246.3422184
##   alpha[9]       283.3462482  2.66369313  277.8474352  281.6341046
##   alpha[10]      219.1562742  2.68745723  213.8511030  217.4226411
##   alpha[11]      258.1861801  2.82267251  252.4833096  256.4499143
##   alpha[12]      228.1907674  2.67821382  223.0705725  226.3842512
##   alpha[13]      242.3797189  2.79698660  237.1069528  240.4358546
##   alpha[14]      268.2588684  2.71732136  262.5688968  266.4373825
##   alpha[15]      242.9239134  2.81060606  237.6852887  241.0575444
##   alpha[16]      245.3760676  2.72130079  240.1230314  243.5545330
##   alpha[17]      232.1473781  2.85795292  226.6159546  230.2257737
##   alpha[18]      240.5631123  2.70221334  235.0794508  238.8687811
##   alpha[19]      253.7390085  2.68243814  248.3902020  252.0113354
##   alpha[20]      241.7204001  2.86416719  236.0308973  239.8677321
##   alpha[21]      248.4956916  2.78723439  243.1073384  246.5444261
##   alpha[22]      225.2258270  2.81762060  219.5195535  223.3754398
##   alpha[23]      228.4431641  2.73461328  222.6051844  226.6854269
##   alpha[24]      245.1508235  2.63501708  240.1351554  243.3596567
##   alpha[25]      234.5263087  2.78128331  228.8581276  232.7482953
##   alpha[26]      253.9806078  2.73452628  248.4030181  252.2410027
##   alpha[27]      254.2860958  2.54808777  249.5777908  252.5757730
##   alpha[28]      242.9664656  2.64060098  237.5545891  241.3097817
##   alpha[29]      217.8840352  2.81530953  212.1043874  216.0266376
##   alpha[30]      241.3974418  2.88412281  235.8968942  239.4616713
##   beta[1]          6.0701699  0.23988837    5.6295599    5.9027579
##   beta[2]          7.0579168  0.25275099    6.5518377    6.8862446
##   beta[3]          6.4848320  0.26529857    5.9752814    6.2994094
##   beta[4]          5.3533972  0.25641266    4.8134621    5.1852477
##   beta[5]          6.5707013  0.24824578    6.0609056    6.4080001
##   beta[6]          6.1633262  0.23565976    5.7099808    6.0048494
##   beta[7]          5.9698319  0.25569841    5.4560169    5.7957586
##   beta[8]          6.4179816  0.25450576    5.9172059    6.2578005
##   beta[9]          7.0513725  0.26580336    6.5452452    6.8709382
##   beta[10]         5.8596212  0.24138057    5.3636823    5.7065060
##   beta[11]         6.8038397  0.25512535    6.2984234    6.6369110
##   beta[12]         6.1222221  0.24219197    5.6385177    5.9814465
##   beta[13]         6.1698860  0.25444626    5.6576177    6.0118561
##   beta[14]         6.6862953  0.25884672    6.1883825    6.5049567
##   beta[15]         5.4226856  0.23724832    4.9625396    5.2628071
##   beta[16]         5.9273245  0.23856428    5.4530161    5.7643915
##   beta[17]         6.2684769  0.24284397    5.7956838    6.1145256
##   beta[18]         5.8443937  0.25438384    5.3403493    5.6773952
##   beta[19]         6.4129012  0.25713723    5.9076234    6.2469126
##   beta[20]         6.0524297  0.23421101    5.6010229    5.8954541
##   beta[21]         6.4014545  0.24194834    5.8873367    6.2503140
##   beta[22]         5.8727762  0.23951242    5.3987479    5.7130497
##   beta[23]         5.7385003  0.25005281    5.2601348    5.5649520
##   beta[24]         5.8906368  0.25106056    5.3839170    5.7282164
##   beta[25]         6.8999612  0.23893841    6.4272533    6.7421982
##   beta[26]         6.5443404  0.24572569    6.0718658    6.3777420
##   beta[27]         5.9074945  0.23014974    5.4374604    5.7517532
##   beta[28]         5.8535090  0.24414947    5.3956719    5.6839514
##   beta[29]         5.6820615  0.24489840    5.2082358    5.5109580
##   beta[30]         6.1387519  0.23648680    5.6925340    5.9719287
##   mu_alpha       242.5178723  2.68611680  237.2267209  240.7470762
##   mu_beta          6.1907515  0.10586261    5.9768566    6.1185673
##   sigmasq_y       37.7777138  5.88761902   27.8452809   33.5348247
##   sigmasq_alpha  214.5946270 61.49038886  124.9580693  169.8228966
##   sigmasq_beta     0.2702446  0.09767191    0.1280761    0.2027839
##   sigma_y          6.1281324  0.47321359    5.2768627    5.7909261
##   sigma_alpha     14.5076174  2.03169713   11.1784636   13.0316111
##   sigma_beta       0.5121260  0.08932866    0.3578772    0.4503154
##   alpha0         106.3213395  3.55335578   99.5836386  103.8685244
##   lp__          -439.0926770  6.64241703 -452.7594430 -443.3948416
##                stats
## parameter                50%          75%        97.5%
##   alpha[1]       239.7953230  241.8226917  245.1152715
##   alpha[2]       247.7401478  249.5451986  252.9416594
##   alpha[3]       252.5223337  254.2400050  257.4981371
##   alpha[4]       232.7265422  234.3602142  237.5521243
##   alpha[5]       231.5837292  233.3976848  236.8976075
##   alpha[6]       249.6935576  251.4262422  255.0387350
##   alpha[7]       228.8267203  230.8046293  234.0596074
##   alpha[8]       248.1138171  250.0857366  253.9541610
##   alpha[9]       283.4004550  285.1599414  288.2149802
##   alpha[10]      219.1309693  220.9094533  224.5082897
##   alpha[11]      258.1539261  260.1009515  263.6489366
##   alpha[12]      228.2154172  229.9946067  233.5722652
##   alpha[13]      242.3646350  244.2949197  247.8227405
##   alpha[14]      268.3734870  270.1208399  273.3090060
##   alpha[15]      242.8941362  244.8184121  248.3579469
##   alpha[16]      245.3857135  247.2367111  250.7617351
##   alpha[17]      232.2593307  234.0816229  237.5371577
##   alpha[18]      240.5258608  242.2339719  246.0926019
##   alpha[19]      253.8017039  255.4674728  259.0434919
##   alpha[20]      241.8444766  243.6287490  247.2336380
##   alpha[21]      248.5045423  250.2830555  254.1742515
##   alpha[22]      225.1548806  227.0718694  230.9600325
##   alpha[23]      228.4623056  230.2020363  233.9244151
##   alpha[24]      245.1735893  246.8749724  250.4368111
##   alpha[25]      234.4709175  236.3795426  240.0058128
##   alpha[26]      254.0052244  255.7341589  259.6545240
##   alpha[27]      254.3213470  255.9263773  259.7263923
##   alpha[28]      242.9339396  244.7512727  247.9803880
##   alpha[29]      217.8812788  219.6377462  223.5253426
##   alpha[30]      241.3461776  243.3199871  247.0525573
##   beta[1]          6.0570743    6.2341566    6.5153920
##   beta[2]          7.0544968    7.2406474    7.5564443
##   beta[3]          6.4859866    6.6679818    6.9935676
##   beta[4]          5.3608791    5.5366693    5.8333519
##   beta[5]          6.5774503    6.7354097    7.0472251
##   beta[6]          6.1621271    6.3159696    6.6412220
##   beta[7]          5.9749113    6.1376668    6.4706856
##   beta[8]          6.4173170    6.5844322    6.9228151
##   beta[9]          7.0489227    7.2211352    7.5893657
##   beta[10]         5.8620173    6.0250047    6.3281135
##   beta[11]         6.8126315    6.9794234    7.3029219
##   beta[12]         6.1214102    6.2774428    6.6014095
##   beta[13]         6.1615392    6.3404230    6.6571696
##   beta[14]         6.6820962    6.8654986    7.1703971
##   beta[15]         5.4274003    5.5803113    5.9101727
##   beta[16]         5.9262474    6.0941318    6.3766716
##   beta[17]         6.2631793    6.4260049    6.7277244
##   beta[18]         5.8406731    6.0096212    6.3422636
##   beta[19]         6.4196736    6.5734420    6.9136626
##   beta[20]         6.0526775    6.2022182    6.5252676
##   beta[21]         6.4127536    6.5604219    6.8741565
##   beta[22]         5.8765874    6.0283719    6.3523970
##   beta[23]         5.7317950    5.9142991    6.2284916
##   beta[24]         5.8945008    6.0516055    6.3634345
##   beta[25]         6.8958102    7.0538883    7.3800807
##   beta[26]         6.5495454    6.7123946    7.0133975
##   beta[27]         5.9087900    6.0611418    6.3470518
##   beta[28]         5.8439617    6.0185156    6.3501240
##   beta[29]         5.6798621    5.8426826    6.1355987
##   beta[30]         6.1407903    6.3061317    6.6028010
##   mu_alpha       242.5758159  244.2478887  247.5292893
##   mu_beta          6.1941616    6.2561782    6.4073906
##   sigmasq_y       37.4916566   41.3133595   51.0177465
##   sigmasq_alpha  203.7956678  247.1597949  361.9708566
##   sigmasq_beta     0.2524830    0.3159583    0.5210212
##   sigma_y          6.1230431    6.4275469    7.1426708
##   sigma_alpha     14.2757020   15.7213166   19.0255287
##   sigma_beta       0.5024766    0.5621017    0.7218180
##   alpha0         106.2670518  108.7476813  113.4326947
##   lp__          -439.0995188 -434.6911579 -426.1944556
## 
## , , chains = chain:4
## 
##                stats
## parameter               mean          sd         2.5%          25%
##   alpha[1]       239.8868051  2.64686717  234.5182418  238.2063914
##   alpha[2]       247.7819390  2.63128207  242.7569624  245.9904534
##   alpha[3]       252.3890607  2.66572282  247.2300563  250.6091994
##   alpha[4]       232.5113827  2.74715931  227.2237380  230.5594109
##   alpha[5]       231.6458566  2.53605373  226.8579344  229.8967219
##   alpha[6]       249.8124634  2.77518481  244.3703152  247.8663916
##   alpha[7]       228.7646347  2.54006449  223.8797812  227.0709479
##   alpha[8]       248.2667698  2.73245767  243.4232575  246.3657291
##   alpha[9]       283.3117646  2.77569442  277.6588566  281.5656057
##   alpha[10]      219.2412076  2.59880496  213.9553585  217.4603954
##   alpha[11]      258.2921712  2.72258137  253.2055014  256.4278565
##   alpha[12]      228.1411716  2.76349840  222.8143739  226.2709706
##   alpha[13]      242.3450226  2.72195205  236.9418525  240.5609815
##   alpha[14]      268.2169669  2.67339446  263.1405667  266.4143105
##   alpha[15]      242.7859829  2.70245770  237.6044974  240.9980768
##   alpha[16]      245.2904552  2.60362460  240.0941027  243.5718959
##   alpha[17]      232.2875888  2.83258217  226.6480463  230.4321124
##   alpha[18]      240.5595852  2.82707994  235.0296508  238.7500433
##   alpha[19]      253.7559447  2.54674784  248.7635450  252.0601418
##   alpha[20]      241.7166269  2.63303689  236.4060387  240.0098012
##   alpha[21]      248.6770200  2.58400847  243.9980593  246.9395452
##   alpha[22]      225.2592710  2.56473381  220.0457234  223.5380858
##   alpha[23]      228.6543331  2.54869403  223.3204606  227.0181786
##   alpha[24]      245.0698896  2.74366813  239.8496764  243.2248450
##   alpha[25]      234.3775581  2.69530248  229.0634876  232.6655846
##   alpha[26]      253.9478189  2.76844900  248.2238886  252.1426363
##   alpha[27]      254.2784326  2.78473953  248.8071095  252.4381398
##   alpha[28]      242.9820091  2.83369821  237.1988045  241.0601964
##   alpha[29]      217.9590093  2.59007205  213.2464454  216.3117548
##   alpha[30]      241.5164540  2.54048044  236.8270761  239.7026616
##   beta[1]          6.0580181  0.24627218    5.5857527    5.9008123
##   beta[2]          7.0507565  0.25620577    6.5189281    6.8812093
##   beta[3]          6.4854125  0.23753143    6.0245502    6.3285150
##   beta[4]          5.3398624  0.24777866    4.8290452    5.1826931
##   beta[5]          6.5732265  0.23704110    6.0949606    6.4124515
##   beta[6]          6.1869150  0.23814487    5.7175528    6.0247618
##   beta[7]          5.9707276  0.23424534    5.5030867    5.8181361
##   beta[8]          6.4185847  0.24347179    5.9416538    6.2553562
##   beta[9]          7.0512614  0.25495670    6.4992986    6.8924669
##   beta[10]         5.8553670  0.24150663    5.3953813    5.6984729
##   beta[11]         6.7932100  0.26659577    6.2873007    6.6201832
##   beta[12]         6.1137323  0.24240055    5.6357968    5.9548262
##   beta[13]         6.1788814  0.25917489    5.6451282    5.9986324
##   beta[14]         6.6926700  0.25305990    6.2142066    6.5259580
##   beta[15]         5.4141822  0.24753669    4.9464824    5.2415233
##   beta[16]         5.9195157  0.25222635    5.4302144    5.7419377
##   beta[17]         6.2773264  0.26090743    5.7459512    6.1085566
##   beta[18]         5.8520853  0.22673299    5.4134053    5.7012770
##   beta[19]         6.3978346  0.23774225    5.9293673    6.2328895
##   beta[20]         6.0502114  0.25365712    5.5867038    5.8746658
##   beta[21]         6.4102521  0.26161467    5.9213165    6.2149224
##   beta[22]         5.8473829  0.24024492    5.3953873    5.6748722
##   beta[23]         5.7499118  0.26039886    5.2631228    5.5749740
##   beta[24]         5.8895225  0.25465043    5.4079135    5.7158520
##   beta[25]         6.9072118  0.25400748    6.4352042    6.7312651
##   beta[26]         6.5508885  0.24454974    6.0663607    6.3979682
##   beta[27]         5.8985094  0.24214781    5.4132212    5.7330473
##   beta[28]         5.8462425  0.23848371    5.3683216    5.6908316
##   beta[29]         5.6739886  0.24597519    5.1870620    5.5182602
##   beta[30]         6.1333733  0.23854329    5.6825510    5.9628392
##   mu_alpha       242.4973392  2.66725037  237.2370232  240.7619799
##   mu_beta          6.1832235  0.11038014    5.9621103    6.1059631
##   sigmasq_y       37.2777193  5.95125437   26.8089044   33.2285181
##   sigmasq_alpha  214.8608747 60.80616835  124.6977843  171.5266041
##   sigmasq_beta     0.2705845  0.10141391    0.1238200    0.2000185
##   sigma_y          6.0867860  0.47852275    5.1777316    5.7644183
##   sigma_alpha     14.5208361  2.00254961   11.1668151   13.0968152
##   sigma_beta       0.5117773  0.09315135    0.3518807    0.4472342
##   alpha0         106.4664232  3.64097096   99.1693945  103.9998035
##   lp__          -438.4282561  7.41528386 -453.3311272 -442.8941137
##                stats
## parameter                50%          75%        97.5%
##   alpha[1]       239.9174372  241.5845400  244.8524258
##   alpha[2]       247.6560853  249.6002015  253.1645567
##   alpha[3]       252.3512595  254.1966649  257.8076243
##   alpha[4]       232.4820886  234.4275281  237.7855349
##   alpha[5]       231.5443774  233.3268598  236.5073489
##   alpha[6]       249.8912492  251.8065450  255.0849111
##   alpha[7]       228.7433729  230.3790644  233.8858232
##   alpha[8]       248.2352212  250.1419205  253.7210835
##   alpha[9]       283.4241544  285.1674916  288.4733256
##   alpha[10]      219.2943985  220.9975449  224.2669334
##   alpha[11]      258.2953720  260.0198678  263.8958646
##   alpha[12]      228.2183191  230.0537923  233.5704259
##   alpha[13]      242.3452957  244.1968484  247.5040710
##   alpha[14]      268.2800079  269.9695423  273.5082773
##   alpha[15]      242.7711752  244.5444550  247.9377072
##   alpha[16]      245.2804796  246.9843015  250.3325518
##   alpha[17]      232.2811229  234.0948592  237.8684299
##   alpha[18]      240.5480460  242.3804493  246.2023281
##   alpha[19]      253.7197835  255.4343062  258.8152436
##   alpha[20]      241.7286483  243.4321528  246.7811340
##   alpha[21]      248.6886099  250.4149197  253.8052086
##   alpha[22]      225.3348428  226.8534969  230.2081357
##   alpha[23]      228.7477778  230.3186177  233.5241400
##   alpha[24]      245.1033207  247.0040683  250.3501972
##   alpha[25]      234.3565560  236.0984075  239.5937414
##   alpha[26]      253.9691438  255.7616071  259.3976438
##   alpha[27]      254.2556889  256.1721080  259.5763990
##   alpha[28]      242.8898041  244.8974316  248.4319710
##   alpha[29]      217.8123672  219.6836418  223.1868211
##   alpha[30]      241.5028561  243.2711098  246.2570559
##   beta[1]          6.0515695    6.2163460    6.5399459
##   beta[2]          7.0570363    7.2293432    7.5131586
##   beta[3]          6.4873006    6.6485816    6.9524299
##   beta[4]          5.3342232    5.4935655    5.8489087
##   beta[5]          6.5720259    6.7352865    7.0303646
##   beta[6]          6.1807716    6.3455350    6.6586552
##   beta[7]          5.9737921    6.1307786    6.4177977
##   beta[8]          6.4112710    6.5856315    6.8964935
##   beta[9]          7.0448015    7.2179397    7.5485346
##   beta[10]         5.8528694    6.0124665    6.3158228
##   beta[11]         6.8005100    6.9565755    7.3361636
##   beta[12]         6.1107271    6.2756795    6.6109746
##   beta[13]         6.1763173    6.3518566    6.6787559
##   beta[14]         6.6836511    6.8541864    7.2174820
##   beta[15]         5.4187277    5.5824149    5.9140072
##   beta[16]         5.9254710    6.1021868    6.3734574
##   beta[17]         6.2892993    6.4540046    6.7626293
##   beta[18]         5.8488151    5.9985340    6.3237274
##   beta[19]         6.4052169    6.5687162    6.8381995
##   beta[20]         6.0458495    6.2312275    6.5549838
##   beta[21]         6.4128112    6.5982596    6.8794340
##   beta[22]         5.8572477    6.0135790    6.3153987
##   beta[23]         5.7372930    5.9283252    6.2819208
##   beta[24]         5.8972562    6.0478470    6.4172307
##   beta[25]         6.9033780    7.0790977    7.4214783
##   beta[26]         6.5489587    6.7149183    7.0411321
##   beta[27]         5.8986770    6.0590043    6.3676260
##   beta[28]         5.8435118    6.0138741    6.3251701
##   beta[29]         5.6840831    5.8287323    6.1522461
##   beta[30]         6.1343568    6.2967186    6.5967581
##   mu_alpha       242.5189834  244.3516290  247.6480769
##   mu_beta          6.1864496    6.2580496    6.3874073
##   sigmasq_y       36.7048170   41.0566691   49.5044001
##   sigmasq_alpha  206.6644391  247.8234694  362.3880683
##   sigmasq_beta     0.2517077    0.3223430    0.5208619
##   sigma_y          6.0584500    6.4075473    7.0359363
##   sigma_alpha     14.3758282   15.7424099   19.0364910
##   sigma_beta       0.5017048    0.5677526    0.7217076
##   alpha0         106.5240975  108.9072703  113.9302470
##   lp__          -438.0246915 -433.3851105 -425.2839243

0.9.4 eight schools example

# stan script
eightschools<-"
// saved as 8schools.stan
data {
  int<lower=0> J;         // number of schools 
  real y[J];              // estimated treatment effects
  real<lower=0> sigma[J]; // standard error of effect estimates 
}
parameters {
  real mu;                // population treatment effect
  real<lower=0> tau;      // standard deviation in treatment effects
  vector[J] eta;          // unscaled deviation from mu by school
}
transformed parameters {
  vector[J] theta = mu + tau * eta;        // school treatment effects
}
model {
  target += normal_lpdf(eta | 0, 1);       // prior log-density
  target += normal_lpdf(y | theta, sigma); // log-likelihood
}"
schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))
fit <- rstan::stan(model_code  = eightschools,
                         data = schools_dat)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
print(fit)
## Inference for Stan model: 432a5acc93a981154b518b22530ee9a6.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu         7.84    0.11 5.01  -2.32   4.58   7.87  11.14  17.95  2145 1.00
## tau        6.63    0.14 5.52   0.34   2.50   5.26   9.24  21.06  1568 1.00
## eta[1]     0.42    0.02 0.95  -1.54  -0.21   0.45   1.04   2.21  3742 1.00
## eta[2]    -0.03    0.01 0.86  -1.74  -0.59  -0.03   0.52   1.65  3510 1.00
## eta[3]    -0.18    0.02 0.91  -1.98  -0.79  -0.20   0.42   1.64  3361 1.00
## eta[4]    -0.02    0.01 0.88  -1.76  -0.58  -0.03   0.53   1.73  3935 1.00
## eta[5]    -0.35    0.02 0.87  -2.04  -0.92  -0.38   0.23   1.37  3198 1.00
## eta[6]    -0.22    0.01 0.89  -1.93  -0.81  -0.22   0.36   1.56  3998 1.00
## eta[7]     0.35    0.01 0.85  -1.37  -0.18   0.37   0.91   2.04  3681 1.00
## eta[8]     0.08    0.02 0.94  -1.77  -0.55   0.09   0.71   1.93  3914 1.00
## theta[1]  11.57    0.16 8.44  -2.06   6.21  10.32  15.83  31.71  2920 1.00
## theta[2]   7.79    0.09 6.26  -4.37   3.88   7.65  11.68  20.55  4550 1.00
## theta[3]   6.16    0.13 7.67 -10.65   2.04   6.59  10.70  21.09  3349 1.00
## theta[4]   7.60    0.10 6.44  -5.40   3.63   7.60  11.62  20.48  3772 1.00
## theta[5]   5.05    0.11 6.27  -8.44   1.21   5.50   9.28  16.06  3451 1.00
## theta[6]   6.16    0.10 6.65  -8.22   2.28   6.61  10.40  18.41  4175 1.00
## theta[7]  10.68    0.12 6.68  -0.78   6.22  10.03  14.55  25.71  3366 1.00
## theta[8]   8.51    0.14 7.91  -7.77   3.98   8.27  12.68  25.18  3414 1.00
## lp__     -39.41    0.07 2.61 -45.24 -40.96 -39.17 -37.59 -34.99  1267 1.01
## 
## Samples were drawn using NUTS(diag_e) at Mon Jan  7 16:18:43 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).
plot(fit)
## 'pars' not specified. Showing first 10 parameters by default.
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

pairs(fit, pars = c("mu", "tau", "lp__"))

la <- extract(fit, permuted = TRUE) # return a list of arrays 
mu <- la$mu 

### return an array of three dimensions: iterations, chains, parameters 
a <- extract(fit, permuted = FALSE) 

### use S3 functions on stanfit objects
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)

0.10 TODO

http://www.stephenpettigrew.com/r/

https://cran.r-project.org/web/packages/apaTables/vignettes/apaTables.html

papaja: Reproducible APA manuscripts with R Markdown

https://crsh.github.io/papaja_man/